Hello R experts,
I have having a difficult time figuring out how to perform and interpret an
ANCOVA of my nested experimental data and would love any suggestions that you
might have.
Here is the deal:
1) I have twelve tanks of fish (1-12), each with a bunch of fish in them
2) I have three treatments (1-3); 4 tanks per treatment. (each tank only has one
treatment applied to it)
3) I sampled multiple fish from each tank (1-3) and would like to nest my tanks
within each treatment (i.e. four tanks nested in treatment 1, four tanks in
treatment 2 and four tanks in treatment 3) in order to account for any
additional variance due to random tank effects within a treatment.
4) The dependent variable in this case is the AREA of an anatomical structure,
which is proportional to body length (the covariate).
5) Here is a simplified example of my data-frame structure and the code to
generate it. I have re-run the following analysis on this example data set for
this email:
treatment<-c(rep(1:3,each=12))
tank<-c(rep(1:12,each=3))
fish<-c(rep(1:3, each=1, times=12))
length<-c(runif(36,15,25))
area<-c(length[1:12]*10+10,length[13:24]*10+20,length[25:36]*10)
fishdata<-data.frame(treatment,tank,fish,length,area);fishdata
# treatment tank fish length area
#1 1 1 1 16.71511 177.1511
#2 1 1 2 21.95281 229.5281
#3 1 1 3 20.39821 213.9821
#4 1 2 1 23.67241 246.7241
#5 1 2 2 17.35663 183.5663
#6 1 2 3 21.61087 226.1087
#7 1 3 1 21.20769 222.0769
#8 1 3 2 20.14224 211.4224
#9 1 3 3 18.50520 195.0520
#10 1 4 1 21.57603 225.7603
#11 1 4 2 20.27112 212.7112
#12 1 4 3 22.78739 237.8739
#13 2 5 1 19.80727 218.0727
#14 2 5 2 15.11602 171.1602
#...and so on...
plot(length,area)
I have tried the following code to run the ANCOVA with the nested design. Area
is the dependent, treatment is a fixed factor, length is a covariate of area,
and tank is nested within treatment.
QUESTION 1. Have I written this code correctly? I a little confused if it should
be tank/treatment of treatment/tank. From what I understand, this part is
testing the equality of the slopes for each treatment group (the interaction
term).
RSslopes<-aov(area~length*treatment+Error(tank/treatment),data=fishdata)
summary(RSslopes)
#Error: tank
# Df Sum Sq Mean Sq
#length 1 753.6 753.6
#
#Error: tank:treatment
# Df Sum Sq Mean Sq
#length 1 4633 4633
#
#Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
#length 1 25282 25282 2254.451 < 2e-16 ***
#treatment 1 327 327 29.192 7.47e-06 ***
#length:treatment 1 5 5 0.461 0.502
#Residuals 30 336 11
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
These results show no significant difference in slopes (P=0.502), so I interpret
that as conforming to the assumption of equal slopes in order to run the actual
ANCOVA and test for a treatment effect.
QUESTION 2. How do I test to see if there is a tank effect at this point,?
I'm not sure which Mean Sq value to test over the residuals, or if I should
wait to do that until the next test?
This next part removes the interaction term to just test the treatment effect
while removing the effect of length. right?
RSancova<-aov(area~length+treatment+Error(tank/treatment),data=fishdata)
summary(RSancova)
#Error: tank
# Df Sum Sq Mean Sq
#length 1 753.6 753.6
#
#Error: tank:treatment
# Df Sum Sq Mean Sq
#length 1 4633 4633
#
#Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
#length 1 25282 25282 2294.34 < 2e-16 ***
#treatment 1 327 327 29.71 5.9e-06 ***
#Residuals 31 342 11
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
So I think this tells me that I do have a treatment effect...but I still would
like to know how to arrange the Mean Sq for the F-test for tank effect.
Now I test to see if removing the interaction term significantly affects the fit
of the model...but this part doesn't work with a nested design!! (it does
work if I do the exact same thing without the nested term.
I get the following:
anova(RSslopes,RSancova)
#Error in UseMethod("anova") :
# no applicable method for 'anova' applied to an object of class
"c('aovlist', 'listof')"
This has something to do with the nested design putting out an aov list, I
think.
Now I really lose it....because I know how to perform a Tukey HSD test to
compare each combination of treatments if I had a regular non-nested ANCOVA, but
when I try it with these tests it can't run...
HSD.test(RSancova,"treatment")
Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors =
stringsAsFactors) :
cannot coerce class 'c("aovlist", "listof")' into
a data.frame
QUESTION 3: any tips at this point to get either of these to work?
I have also tried another method that I found on the help list archive:
require(MASS) ## for oats data set
require(nlme)## for lme()
require(multcomp)## for multiple comparison stuff
Aov.mod <- aov(area~ length.mm+ treatment + Error(tank/treatment), data =
d.sag.right)
Lme.mod <- lme(area ~ length.mm+ treatment, random = ~1 | tank/treatment,
data = d.sag.right)
summary(Aov.mod)
#Error: tank
# Df Sum Sq Mean Sq
#length.mm 1 8181085 8181085
#
#Error: tank:treatment
# Df Sum Sq Mean Sq
#length.mm 1 1822796 1822796
#treatment 1 11304871 11304871
#
#Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
#length.mm 1 66978274 66978274 161.96 < 2e-16 ***
#treatment 2 8871292 4435646 10.73 0.000106 ***
#Residuals 59 24399419 413549
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(Lme.mod)
# numDF denDF F-value p-value
#(Intercept) 1 54 4207.021 <.0001
#length.mm 1 54 183.772 <.0001
#treatment 2 8 9.779 0.0071
summary(Lme.mod)
#Linear mixed-effects model fit by REML
# Data: d.sag.right
# AIC BIC logLik
# 1011.271 1026.16 -498.6353
#
#Random effects:
# Formula: ~1 | tank
# (Intercept)
#StdDev: 193.3108
#
# Formula: ~1 | treatment %in% tank
# (Intercept) Residual
#StdDev: 193.3063 635.9044
#
#Fixed effects: area ~ length.mm + treatment
# Value Std.Error DF t-value p-value
#(Intercept) -3005.7828 741.2499 54 -4.055019 0.0002
#length.mm 523.9419 37.1477 54 14.104286 0.0000
#treatment1000 588.6912 288.0330 8 2.043832 0.0752
#treatment2000 1283.1145 292.3259 8 4.389329 0.0023
# Correlation:
# (Intr) lngth. tr1000
#length.mm -0.956
#treatment1000 -0.246 0.025
#treatment2000 -0.384 0.173 0.567
#
#Standardized Within-Group Residuals:
# Min Q1 Med Q3 Max
#-2.487003437 -0.594427676 0.004423702 0.445992447 2.463725291
#
#Number of Observations: 66
#Number of Groups:
# tank treatment %in% tank
# 11 11
summary(glht(Lme.mod, linfct=mcp("treatment"="Tukey")))
*** caught bus error ***
address 0x0, cause 'non-existent physical address'
etc. R crash notification...
This part always crashes R...
QUESTION 4. WTF? why does it crash R?
QUESTION 5. How do I interpret the anova and summary of the Lme.mod?
QUESTION 6. Why does the Number of Groups at the bottom say 11? and say
treatment %in% tank...that makes me think that I didn't nest it properly.
Well...I think that is it...I'm confused. please help!!
--Sean
[[alternative HTML version deleted]]
