假设您有这样的数据:
Y <- rnorm(100)
x <- rnorm(100)
z <- rnorm(100)
data <- data.frame(Y, x, z)
为此，我们可以拟合线性模型:
fit <- lm(Y ~ x + z, data)
文章指出，回归估计是通过求解一组正态方程来进行的:
或矩阵形式:
有了lm个对象，我如何提取d和S矩阵呢？
假设您有这样的数据:
Y <- rnorm(100)
x <- rnorm(100)
z <- rnorm(100)
data <- data.frame(Y, x, z)
为此，我们可以拟合线性模型:
fit <- lm(Y ~ x + z, data)
文章指出，回归估计是通过求解一组正态方程来进行的:
或矩阵形式:
有了lm个对象，我如何提取d和S矩阵呢？
我们可以从fm得到X和Y，如图所示.然后用crossprod计算XtX和XtY，然后用solve计算系数.
fm <- lm(mpg ~ cyl + carb, mtcars)
X <- model.matrix(fm)
Y <- model.response(model.frame(fm))
# or
Y <- fitted(fm) + resid(fm)
XtX <- crossprod(X); XtX
## (Intercept) cyl carb
## (Intercept) 32 198 90
## cyl 198 1324 604
## carb 90 604 334
XtY <- crossprod(X,Y); XtY
## [,1]
## (Intercept) 642.9
## cyl 3693.6
## carb 1641.9
c(solve(XtX, XtY))
## [1] 37.812739 -2.625023 -0.526146
coef(fm)
## (Intercept) cyl carb
## 37.812739 -2.625023 -0.526146