在R中重新创建SAS混合模型输出（包括F测试）

Posted 2021-04-20

我最近在SAS学习了ANOVA课程，并在R中重写了我的代码。到目前为止，将随机效果（和混合效果）模型从SAS转换为R已经躲过了我。我从R得到的输出与SAS非常不同：SS和F值不同，我无法对随机效应进行F检验。我能得到的最接近的是Chi-sq，使用rand（）。所以也许我在R中做错了。

以下是SAS代码和输出，然后是我在R中的尝试。

dataset download, <1kb

*Two-Way ANOVA, with one random effect and interaction term;
*import dataset as "pesticide";
proc glm data=pesticide;
    class locations chemicals;
    model numberkilled = locations chemicals locations*chemicals / solution;
    random locations locations*chemicals / test;
run; quit;

以下是尝试的R代码。

#data step
pesticide <- read.csv("ex17-10.txt")
colnames(pesticide) <- c("location", "chemical", "number_killed")
pesticide$location <- as.factor(pesticide$location)
pesticide$chemical <- as.factor(pesticide$chemical)

#ANOVA
library(lmerTest); library(car)
model <- lmer(number_killed ~ chemical + (1|location) + (1|chemical:location), data=pesticide)
Anova(model, type=3, test="F")

输出是下一个。随机效应和相互作用项（也是随机的）没有F检验，SS和F值与SAS不同。

Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: number_killed
                  F Df Df.res    Pr(>F)    
(Intercept) 587.069  1     16 4.879e-14 ***
chemical     48.108  3     12 5.800e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

总之，我不知道如何在R中正确地做混合效果模型。固定效果模型都可以。

答案

您可以按照here给出的建议重现SAS代码的输出：

fit <- lm(number_killed ~ location * chemical, data=pesticide)
results <- anova(fit)
Df <- results$Df
SumSq <- results$"Sum Sq"
MeanSq <- results$"Mean Sq"
Fvalue <- results$"F value"
Pvalue <- results$"Pr(>F)"
Error.Term <- MeanSq[3]
df.error <- Df[3]

Fvalue[1] <- MeanSq[1]/Error.Term
Pvalue[1] <- 1 - pf(Fvalue[1], Df[1], df.error)

Fvalue[2] <- MeanSq[2]/Error.Term
Pvalue[2] <- 1 - pf(Fvalue[2], Df[2], df.error)

Ftable <- cbind(Df, SumSq, MeanSq, Fvalue, Pvalue)
rownames(Ftable) <- c("Locations", "Chemicals", "Locations:Chemicals", "Residuals")
print(Ftable)

#                     Df    SumSq    MeanSq     Fvalue       Pvalue
# Locations            4   3.8115  0.952875  0.7076461 6.020037e-01
# Chemicals            3 180.1327 60.044250 44.5914534 8.797523e-07
# Locations:Chemicals 12  16.1585  1.346542  3.8889290 3.652306e-03
# Residuals           20   6.9250  0.346250         NA           NA