模型中嵌套效应时JMP与R之间的平方和差异

js4nwp54 于 2023-03-10 发布在其他

关注(0)|答案(2)|浏览(172)

我尝试将JMP脚本转换为R。在我的分析中，我必须执行以下线性模型：Y ~ G + C + G*C + C%in%L（在jmp中编写的等效模型将是G+ C + G * C + L[C]）。
我知道JMP给出的是类型3 SSQ，我也知道如何使用适当的对比在R中获得等价。但是，我认为由于我有这种嵌套效应，有一个项在R和JMP中仍然有不同的SSQ，我不知道为什么以及如何获得相同的值。下面是一些模拟来说明这个问题：

# R code
# library('data.table')
set.seed(18)
TEST=data.table(
 L= c( rep( LETTERS[seq( from = 1, to = 3 )],3), rep( LETTERS[seq( from = 4, to = 6 )],3),rep( LETTERS[seq( from = 7, to = 9 )],3)),
 G= rep(sort(rep(LETTERS[seq( from = 1, to = 3 )],3)),3),
 C=sort(rep(LETTERS[seq( from = 10, to = 12)],9)),
 Y= rnorm(27)
 )
 set.seed(18) 
# part add to have the data not orthogonal
TEST_add=data.table(
 L= sort(rep( LETTERS[seq( from = 10, to = 14 )],2)),
G= rep(LETTERS[seq( from = 1, to = 2 )],5),
C = rep("L",10),
 Y= rnorm(10)
)
TEST=rbind(TEST_add,TEST)
 
 TEST$L = as.factor(TEST$L)
 TEST$G = as.factor(TEST$G)
 TEST$C = as.factor(TEST$C)  
   
model <- lm(Y~ G  +C+ C*G +  C%in% L, data = TEST ,contrasts=list(G="contr.sum",C="contr.sum", L="contr.sum")) 
drop1(model,.~.,test="F",all.cols = FALSE)

这将在R中给予此表：

Df Sum of Sq    RSS    AIC F    value Pr(>F)
<none>                24.196   24.285               
G       2    0.2945   24.491   20.733  0.1035 0.9023
C       2    0.4707   24.667   20.998  0.1653 0.8490
G:C     4    5.9115   30.108   24.373  1.0383 0.4163
C:L    11   21.7349   45.931   26.000  1.3882 0.2631

但JMP将返回此表

Source  DF  Sum of Squares
C       2   0.876575
G       2   0.294510
L[C]    11  21.734947
G*C     4   5.911491

一切都很好，除了“C”：在JMP中SSQ为0.877，在R中为0.47
我试过改变对比，但问题似乎不在这里。当我只有固定效应或简单交叉效应时，我有完全相同的数字。我猜嵌套泰尔梅在两个程序中没有以相同的方式考虑？

来源：https://stackoverflow.com/questions/75594861/differences-in-sum-of-squares-between-jmp-and-r-when-nested-effects-in-the-model

2条答案

按热度按时间

htzpubme1#

鉴于您指出sasLM可以重现SAS III型平方和，因此我们完全修改了此答案，以便进行良好的比较。
III型平方和通常表示为不含感兴趣项的模型的残差平方和与完整模型的残差平方和之差;然而，SAS类型III的工作方式有些不同。
代替将模型矩阵列划分成两个集合，它将模型矩阵列划分成3个集合：

X 0是不包含感兴趣的主效应的列，
X1是感兴趣的主效应的列，
X2是表示包含该主效应的交互作用的列

为了获得平方和，它将由X 0和X2 s组成的模型与完整模型进行比较，其中X2 s等于以下X 0、X1和X2项。

X2s <- X2 %*% t(X2) %*% resid(lm(X2 ~ X0 + X1 + 0, TEST))

就问题中的示例而言，假设我们对主效应C平方和感兴趣，首先创建一个模型矩阵X，没有对比度。查看colnames(X)，我们可以看到X0，X1和X2由以下代码中指示的列定义。然后定义X2s，并取包含X0和X2s的模型的离差与完整模型的离差之差。请注意，deviance是残差平方和（对于线性模型）。有关所有这些计算背后的线性代数的解释，请参见LaMotte (2018)。

fo <- Y ~ G + C + C * G + C/L
X <- model.matrix(fo, TEST, contrasts.arg = 
  lapply(TEST[, c("G", "C", "L")], contrasts, contrasts = FALSE))

X0 <- X[, 1:4]  # cols not involving C
X1 <- X[, 5:7]  # C main effect
X2 <- X[, 8:58] # columns containing C other than X1

X2s <- X2 %*% t(X2) %*% resid(lm(X2 ~ X0 + X1 + 0, TEST))
deviance(lm(Y ~ X0 + X2s + 0, TEST)) - deviance(lm(Y ~ X + 0, TEST))
## [1] 0.8765751  <--- same as sasLM sum of squares for C

赞(0）回复(0）举报 2023-03-10

pbpqsu0x2#

我的同事确实为我找到了一个解决方案，我在这里分享它，以防它对其他人有用。R包'sasLM'允许使用JMP（或SAS）中的相同值。

library(sasLM)
aov3(Formula = Y~ G  +C+ C*G +  C/L, Data=TEST, BETA=FALSE, Resid=FALSE)

是否返回

Response : Y
                Df Sum Sq Mean Sq F value Pr(>F)
MODEL           19 30.508 1.60568  1.1281 0.4041
 G               2  0.295 0.14725  0.1035 0.9023
 C               2  0.877 0.43829  0.3079 0.7390
 G:C             4  5.911 1.47787  1.0383 0.4163
 C:L            11 21.735 1.97590  1.3882 0.2631
RESIDUALS       17 24.196 1.42331               
CORRECTED TOTAL 36 54.704

完全吻合！
尽管如此，我还是不知道为什么结果会不同，我对答案很感兴趣！

赞(0）回复(0）举报 2023-03-10

我来回答

模型中嵌套效应时JMP与R之间的平方和差异

2条答案

相关问题

热门标签

最新问答