CODEWITHSUNDEEP
rsasanovamixed-modelsrandom-effects
dkae
dkae

Reputation: 313

Recreating SAS mixed model output (including F tests) in R

I recently took an ANOVA class in SAS, and am rewriting my code in R. Thus far, translating random effect (and mixed effect) models from SAS to R has eluded me. The output I get from R is very different from SAS: the SS and F value are different, and I can't get F tests for the random effects. The closest I've been able to get are Chi-sq, using rand(). So perhaps I'm doing it all wrong in R.

The following is the SAS code and output, and then the attempt I made in 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;

SAS output

The following is the attempted R code.

#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")

Output is next. There are no F tests for the random effect and interaction term (which is also random), and the SS and F value are different from 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

In summary, I'm don't know how to properly do mixed effect models in R. Fixed effects models are all ok.

Upvotes: 2

Views: 752

Answers (1)

Marco Sandri
Marco Sandri

Reputation: 24272

You can reproduce the output of the SAS code following the suggestions given here:

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

Upvotes: 1

Related Questions