Re: [R-sig-eco] ANCOVA with random effects for slope and intercept
peterhouk1 . peterhouk@... 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 BIClogLik 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 Residual3.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.71008 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 BIClogLik 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 Residual3.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.71008 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
Re: [R-sig-eco] ANCOVA with random effects for slope and intercept
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 . peterh...@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.61643 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.0954 -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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology [[alternative HTML version deleted]] ___ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Re: [R-sig-eco] ANCOVA with random effects for slope and intercept
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? Mlme1-lme(response ~ predictor, random = ~1 + predictor | group_factor, data=mydata) Linear mixed-effects model fit by REML Data: mydata AIC BIClogLik 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 Residual3.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.71008 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 BIClogLik 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 Residual3.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.71008 2.136060167 Number of Observations: 72 Number of Groups: 11 On Tue, Nov 4, 2014 at 10:41 PM, Mollie Brooks mbro...@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 . peterh...@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
