require(nlme) summary(RatPupWeight) # weight sex Litter Lsize Treatment # Min. :3.680 Male :171 7 : 18 Min. : 2.00 Control:131 # 1st Qu.:5.650 Female:151 9 : 17 1st Qu.:12.00 Low :126 # Median :6.055 8 : 17 Median :14.00 High : 65 # Mean :6.081 11 : 16 Mean :13.33 # 3rd Qu.:6.397 20 : 16 3rd Qu.:16.00 # Max. :8.330 14 : 15 Max. :18.00 (Other):223 plot(RatPupWeight) ## default dot plot plot.design(RatPupWeight[,-4]) ## shows spreads of each factor variable require(lattice) dotplot(Litter ~ weight,group=interaction(Treatment,sex), data =RatPupWeight, autoKey=T ) with(RatPupWeight, interaction.plot( Treatment,sex,weight)) with(RatPupWeight, tapply( weight, interaction(Treatment,sex),summary)) #$Control.Male # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.570 6.060 6.410 6.471 6.940 8.330 #$Low.Male # Min. 1st Qu. Median Mean 3rd Qu. Max. # 5.250 5.750 6.000 6.025 6.240 7.130 #$High.Male # Min. 1st Qu. Median Mean 3rd Qu. Max. # 5.010 5.360 5.690 5.918 6.290 7.700 #$Control.Female # Min. 1st Qu. Median Mean 3rd Qu. Max. # 3.680 5.822 6.175 6.116 6.448 7.570 #$Low.Female # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.750 5.600 5.840 5.838 6.110 7.730 #$High.Female # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.480 5.415 5.760 5.852 6.240 7.680 round(with(RatPupWeight, tapply( weight, interaction(Treatment,sex),sd)),3) # Control.Male Low.Male High.Male Control.Female Low.Female # 0.754 0.380 0.691 0.685 0.450 # High.Female # 0.600 rat.fit1 <- lme(weight ~ Treatment * sex +Lsize, data = RatPupWeight, random = ~1|Litter) rat.fit0 <- lm(weight ~ Treatment * sex + Lsize, data = RatPupWeight) anova(rat.fit1,rat.fit0) # Model df AIC BIC logLik Test L.Ratio p-value #rat.fit1 1 9 421.3015 455.0747 -201.6508 #rat.fit0 2 8 508.7072 538.7277 -246.3536 1 vs 2 89.40562 <.0001 ## or look at the intervals command output to see how close the ## variance estimate for Litter is to zero:: intervals(rat.fit1) #Approximate 95% confidence intervals # # Fixed effects: # lower est. upper #(Intercept) 7.3994490 7.86564125 8.33183350 #Treatment.L -0.9208598 -0.64067890 -0.36049802 #Treatment.Q -0.2388612 0.01144004 0.26174130 #sexFemale -0.4479878 -0.34805818 -0.24812854 #Lsize -0.1671764 -0.12838210 -0.08958782 #Treatment.L:sexFemale -0.1076943 0.07567678 0.25904790 #Treatment.Q:sexFemale -0.1866598 -0.02478407 0.13709169 # Random Effects: # Level: Litter # lower est. upper #sd((Intercept)) 0.22266 0.31067 0.43346 #### Not too close to zero ##### # Within-group standard error: # lower est. upper #0.3728742 0.4043370 0.4384546 plot(rat.fit1, col=as.numeric(RatPupWeight$Treatment)) rat.fit2 <- update(rat.fit1, weights=varPower()) anova(rat.fit1,rat.fit2) # Model df AIC BIC logLik Test L.Ratio p-value #rat.fit1 1 9 421.3015 455.0747 -201.6508 #rat.fit2 2 10 362.6271 400.1529 -171.3136 1 vs 2 60.67441 <.0001 plot(rat.fit2, col=as.numeric(RatPupWeight$Treatment)) rat.fit3 <- update(rat.fit1, weights=varIdent( form =~1|Treatment)) rat.fit4 <- update(rat.fit1, weights=varIdent( form =~1|factor(Treatment=="Control"))) anova(rat.fit1,rat.fit4, rat.fit3) # Model df AIC BIC logLik Test L.Ratio p-value #rat.fit1 1 9 421.3015 455.0747 -201.6508 #rat.fit4 2 10 383.2780 420.8037 -181.6390 1 vs 2 40.02358 <.0001 #rat.fit3 3 11 384.0819 425.3602 -181.0410 2 vs 3 1.19605 0.2741 AIC(rat.fit2,rat.fit3, rat.fit4) # df AIC #rat.fit2 10 362.6271 #rat.fit3 11 384.0819 #rat.fit4 10 383.2780 anova(rat.fit2) # numDF denDF F-value p-value #(Intercept) 1 292 8603.901 <.0001 #Treatment 2 23 4.456 0.0231 #sex 1 292 57.027 <.0001 #Lsize 1 23 34.075 <.0001 #Treatment:sex 2 292 0.000 0.9996 ### Note: order changes these p-values. ### First put in variables an expert would say should be in the model. ### Treatment should go in after the other effects have been accounted for. anova(update(rat.fit2, .~Lsize+sex*Treatment)) # numDF denDF F-value p-value #(Intercept) 1 292 8603.901 <.0001 #Lsize 1 23 17.773 0.0003 #sex 1 292 59.994 <.0001 #Treatment 2 23 11.123 0.0004 #sex:Treatment 2 292 0.000 0.9996 rat.fit5 <- update(rat.fit2, fixed=.~Lsize + Treatment + sex) ## wrong comparison: anova(rat.fit5, rat.fit2) # Model df AIC BIC logLik Test L.Ratio p-value #rat.fit5 1 8 351.6799 381.7511 -167.8399 #rat.fit2 2 10 362.6271 400.1529 -171.3136 1 vs 2 6.947232 0.031 #Warning message: #Fitted objects with different fixed effects. REML comparisons are not meaningful. # in: anova.lme(rat.fit5, rat.fit2) ## right: compare fixed effects with ML estimates anova(update(rat.fit5, method="ML"), update(rat.fit2, method="ML")) # Model df AIC BIC logLik Test #update(rat.fit5, method = "ML") 1 8 332.0872 362.2836 -158.0436 #update(rat.fit2, method = "ML") 2 10 336.0461 373.7916 -158.0230 1 vs 2 # L.Ratio p-value #update(rat.fit5, method = "ML") #update(rat.fit2, method = "ML") 0.041152 0.9796 ## check residuals for normality qqnorm(resid(rat.fit5,typ="p")) qqline(resid(rat.fit5,typ="p")) ## can also look at the random effects which should be normally distributed: qqnorm(ranef(rat.fit5)[,1]) qqline(ranef(rat.fit5)[,1])