1

我有这个数据

y x1 x2 pre
1  16  1  1  14
2  15  1  1  13
3  14  1  2  14
4  13  1  2  13
5  12  2  1  12
6  11  2  1  12
7  11  2  2  13
8  13  2  2  13
9  10  3  1  10
10 11  3  1  11
11 11  3  2  11
12  9  3  2  10

我安装了以下模型

lm(y ~ x1 + x2 + x1*x2)

我的设计矩阵是

  [,1] [,2] [,3] [,4] [,5] [,6] [,7]
 [1,]    1   14    1    0    1    1    0
 [2,]    1   13    1    0    1    1    0
 [3,]    1   14    1    0    0    0    0
 [4,]    1   13    1    0    0    0    0
 [5,]    1   12    0    1    1    0    1
 [6,]    1   12    0    1    1    0    1
 [7,]    1   13    0    1    0    0    0
 [8,]    1   13    0    1    0    0    0
 [9,]    1   10    0    0    1    0    0
[10,]    1   11    0    0    1    0    0
[11,]    1   11    0    0    0    0    0
[12,]    1   10    0    0    0    0    0

我正在尝试使用此设计来重现下表：

Source DF Squares Mean Square F Value Pr > F
Model 6 44.79166667 7.46527778 12.98 0.0064
Error 5 2.87500000 0.57500000
Corrected Total 11 47.66666667

Source DF Type III SS Mean Square F Value Pr > F
pre 1 3.12500000 3.12500000 5.43 0.0671
x1 2 4.58064516 2.29032258 3.98 0.0923
x2 1 3.01785714 3.01785714 5.25 0.0706
x1*x2 2 1.25000000 0.62500000 1.09 0.4055

第一部分很好

XtX <- t(x) %*% x
XtXinv <- solve(XtX)
betahat <- XtXinv %*% t(x) %*% y

H <- x %*% XtXinv %*% t(x) 
IH <- (diag(1,12) - H)
yhat <- H %*% y 
e <- IH %*% y
ybar <- mean(y)

MSS <- t(betahat) %*% t(x) %*% y - length(y)*(ybar^2) 
ESS <- t(e) %*% e 
TSS <- MSS + ESS 

dfM <- sum(diag(H)) - 1 
dfE <- sum(diag(IH)) 
dfT <- dfM + dfE 

MSM <- MSS/dfM 
MSE <- ESS/dfE 

Ftest <- MSM / MSE
pr <- 1 - pf(Ftest, dfM, dfE)

'pre' 的对比度系数矩阵似乎是正确的。

L <- matrix(c(0,1,0,0,0,0,0), 1, 7, byrow=T)
Lb <- L %*% betahat 
LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4) 
SSpre <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb) 
MSpre <- SSpre / 1 
Fpre <- MSpre / MSE 
PRpre <- 1 - pf(Fpre, 1, 12-7)

但我不明白如何定义 x1、x2 和 x1*x2 的对比度系数矩阵。我的其余代码有什么问题？下面是我认为应该如何计算 x1 的示例

L <- matrix(c(0,0,1,1,0,0,0), 1, 7, byrow=T)
Lb <- L %*% betahat 
LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4) 
SSX1 <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb) 
MSX1 <- SSX1 / 1 
FX1 <- MSX1 / MSE 
PRX1 <- 1 - pf(FX1, 1, 12-7)

谢谢！

4

0 回答 0