[R-sig-eco] ANCOVA with random effects for slope and intercept
Ben Bolker
bbolker at gmail.com
Thu Nov 6 03:15:46 CET 2014
peterhouk1 . <peterhouk at ...> writes:
>
> Greetings Mollie -
>
> Sure, the first general approach without explicitly telling R of my
> grouping factor (not sure if that makes a difference, but second example
> below does this). My best guess is that the full model has significance,
> but the random effects model does not. However, simple partial
> correlations show that within many grouping factors
> the relationship holds,
> but not in all. If this were the case, why would our Corr = 0?
There are a few things going on here.
* Plotting your data (see below), it looks pretty clear that
you've already mean-corrected it (all groups have mean value
of zero at the center point of the data) -- you probably want
to take this into account in your model.
* Even with this correction, you still get an estimate of zero
for the variance in slopes among groups. This is not terribly
surprising for small-to-moderate sized, fairly
noisy (large residual variance) data sets. This does indeed
mean that you will get essentially the same answer from a
plain old linear model.
* I did the stuff below with lmer, but I'd expect similar answers
from lme . I used lmerTest for easy access to denominator df
approximates.
* The overall effect is not super-strong (R^2 ~ 0.12), but
very clearly statistically significant (depending on how you
do the calculation, p-value ranges from 0.001 to 0.009 ...)
## get data, draw pictures
mydata <- read.table(header=TRUE,"houk.txt")
library("ggplot2"); theme_set(theme_bw())
## make grouping variable into a factor; not totally necessary
## but doesn't hurt/probably a good idea
mydata <- transform(mydata,group_factor=factor(group_factor))
meanpred <- mean(mydata$predictor)
ggplot(mydata,aes(predictor,response,colour=group_factor))+geom_point()+
geom_smooth(method="lm",se=FALSE)+
geom_vline(xintercept=meanpred)
## centered version of predictor
mydata <- transform(mydata,cpred=predictor-meanpred)
library("lme4")
## original random-slopes model
(m1 <- lmer(response~predictor+(predictor|group_factor), data=mydata))
## model with centered predictor; suppress among-group variation in
## intercept
(m2 <- lmer(response~cpred+(cpred-1|group_factor), data=mydata))
summary(m2B <- lm(response~cpred, data=mydata))
## get p-value, using common-sense/classical value of 10 for denominator df
(cc1 <- coef(summary(m2)))
tval <- cc1["cpred","t value"]
pt(abs(tval),df=10,lower.tail=FALSE)
## refit model with lmerTest for convenient access to df estimates
detach("package:lme4")
library("lmerTest")
(m3 <- lmer(response~cpred+(cpred-1|group_factor), data=mydata))
anova(m3,ddf="Kenward-Roger") ## gives denDF=8.67, p-value=0.0097
anova(m3,ddf="Satterthwaite") ## gives denDF=69.9 (?), p-value=0.001
## essentially equivalent to lm() because denDF are estimated to be large
detach("package:lmerTest")
>
> Mlme1<-lme(response ~ predictor,
> random = ~1 + predictor | group_factor, data=mydata)
>
> Linear mixed-effects model fit by REML
> Data: mydata
> AIC BIC logLik
> 74.80524 88.29622 -31.40262
>
> Random effects:
> Formula: ~1 + predictor | group_factor
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 1.106112e-05 (Intr)
> predictor 1.577405e-10 0
> Residual 3.492176e-01
>
> Fixed effects: response ~ predictor
> Value Std.Error DF t-value p-value
> (Intercept) 0.29852308 0.09774997 60 3.053945 0.0034
> predictor -0.03258404 0.00970759 60 -3.356554 0.0014
> Correlation:
> (Intr)
> predictor -0.907
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -2.560560620 -0.688713759 -0.008759271 0.710084444 2.136060167
>
> Number of Observations: 72
> Number of Groups: 11
>
> Second example where I explicitly tell R of the grouping factor:
>
> mydata$fgroup_factor <- factor(mydata$group_factor)
>
> Mlme1<-lme(response ~ predictor,
> random = ~1 + predictor | fgroup_factor, data=mydata)
>
> Linear mixed-effects model fit by REML
> Data: mydata
> AIC BIC logLik
> 74.80524 88.29622 -31.40262
>
> Random effects:
> Formula: ~1 + predictor | fgroup_factor
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 1.106112e-05 (Intr)
> predictor 1.577405e-10 0
> Residual 3.492176e-01
>
> Fixed effects: response ~ predictor
> Value Std.Error DF t-value p-value
> (Intercept) 0.29852308 0.09774997 60 3.053945 0.0034
> predictor -0.03258404 0.00970759 60 -3.356554 0.0014
> Correlation:
> (Intr)
> predictor -0.907
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -2.560560620 -0.688713759 -0.008759271 0.710084444 2.136060167
>
> Number of Observations: 72
> Number of Groups: 11
>
> On Tue, Nov 4, 2014 at 10:41 PM, Mollie Brooks <mbrooks <at> ufl.edu> wrote:
>
> > Hi Dr. Houk,
> >
> > You say you get the same result from the lmer model as a linear model.
> > Can you show us the summary of both models so that we might help you
> > interpret it?
> >
> > Thanks,
> > Mollie
> > ------------------------
> > Mollie Brooks, PhD
> > Postdoctoral Researcher, Population Ecology Research Group
> > Institute of Evolutionary Biology & Environmental Studies, University of
> > Zürich
> > http://www.popecol.org/team/mollie-brooks/
> >
> >
> > On 4Nov 2014, at 13:30, peterhouk1 . <peterhouk <at> gmail.com> wrote:
> >
> > Greetings -
> >
> > Looking for advice and insight into ANCOVA models that allow for random
> > slope and intercept effects. I have been using lme and lmer, but I can't
> > seem to figure out where I've gone wrong. I have a grouping factor,
> > predictor, and response. The response~predictor relationship should be
> > nested within grouping factor, with a desire to allow for random effects of
> > the slope and y-intercept. Data and my code are found below, greatly
> > appreciate insight.
> >
> > I get the same response from both approaches:
> >
> > Mlme1<-lme(response~predictor,
> > random = list(group_factor = ~1 | predictor), data=mydata)
> >
> > Second approach
> >
> > Mlmer1<-lmer(response~predictor + (1 + predictor|group_factor),
> > data=mydata)
> >
> > Taking both approaches, I get the same results as a simple linear model
> > that does not account for nesting within the group factor.
> >
> > mydata
> >
> > group_factor predictor response 1 13.75584744 -0.257259794 1 3.059971584
> > 0.703113472 1 14.00447296 -0.260892287 1 6.705509001 -0.33269593 1
> > 6.592728067 0.053814446 1 9.211047122 -0.002776485 2 12.50497696
> > 0.22727311 2 7.059077939 0.18586719 2 6.617249805 -0.022951714 2
> > 1.559557719 0.833397702 2 12.1962121 -0.57159955 2 11.29647955
> > -0.963754818 2 15.54334219 0.489700014 3 8.93626518 -0.421471799 3
> > 1.681675438 0.383260174 3 13.43826892 -0.330882653 3 10.8971089
> > 0.47675377
> > 3 7.600869443 -0.227033926 3 15.004137 0.104257183 4 12.61327214
> > 0.131460302 4 3.788474342 0.079758467 4 14.37299098 0.294254076 4
> > 3.564225024 -0.006581881 4 13.63301652 -0.498890965 5 4.36777119
> > 0.90215334 5 5.150130473 0.098069832 5 7.875920526 0.468166528 5
> > 13.91344257 -0.291551635 5 10.1061938 -0.35162982 5 13.32810817
> > 0.133439251 5 15.64845612 -0.439418202 5 1.857959976 -0.468857357 5
> > 11.31495202 -0.050371938 6 7.851162116 0.126588358 6 5.285251391
> > 0.212699384 6 15.82353883 -0.202005195 6 11.90209318 0.34412633 6
> > 5.547563146 -0.446233668 6 5.645270991 -0.303913602 6 8.668938138
> > 0.268738394 7 15.66063395 -0.194882955 7 11.11830972 0.196309336 7
> > 3.174056086 0.188199892 7 12.39006821 0.222261643 7 3.034014836
> > 0.039201594 7 13.13551529 -0.451089511 8 9.990570964 -0.128411654 8
> > 2.382739273 0.145897042 8 4.870002628 0.875223363 8 5.99541207
> > 0.155610264
> > 8 16.56247927 -0.461697164 8 14.58286908 -0.514244888 8 12.72137648
> > -0.072376963 9 7.436676598 -0.306969441 9 5.196390457 -0.222063871 9
> > 7.213670545 -0.411344562 9 9.243636562 0.095582355 9 7.029863529
> > 0.31586562 9 6.074354561 0.459471079 9 8.282168564 0.632851023 9
> > 13.85105359 -0.298326302 9 5.972744894 -0.180974348 9 14.09541111
> > -0.084091553 10 5.304552826 0.265618935 10 4.089248332 0.270559629 10
> > 7.482418571 0.049764287 10 16.44388769 0.008163962 10 13.78009474
> > -0.594106812 11 10.34249094 -0.206619307 11 4.491492572 -0.207263365 11
> > 7.350436798 0.083703414 11 8.38636562 0.330179258
> >
> > --
> >
> > Peter Houk, PhD
> > Assistant Professor
> > University of Guam Marine Laboratory
> > http://www.guammarinelab.com/peterhouk.html
> > www.pacmares.com
> > www.micronesianfishing.com
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-ecology mailing list
> > R-sig-ecology <at> r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
> >
> >
> >
>
More information about the R-sig-ecology
mailing list