[R] AR1 covariance structure for lme object and R/SAS differences in model output
Dear R users,
We are working on a data set in which we have measured repeatedly a
physiological response variable (y)
every 20 min for 12 h (time variable; 'x') in subjects ('id') beloning to
one of five groups ('group'; 'A' to 'E'). Data are located at:
https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0
We are interested to model if the response in y differences with time (i.e.
'x') for the two groups. Thus:
require(nlme)
m1<-lme(y~group*x+group*I(x^2),random=~x|id,data=data.df,na.action=na.omit)
But because data are collected repeatedly over short time intervals for each
subject, it seemed prudent to consider an autoregressive covariance
structure. Thus:
m2<-update(m1,~.,corr=corCAR1(form=~x|id))
AIC values indicate the latter (i.e. m2) as more appropriate:
anova(m1,m2)
# Model df AIC BIC logLikTest L.Ratio
p-value
#m1 1 19 2155.996 2260.767 -1058.9981
#m2 2 20 2021.944 2132.229 -990.9718 1 vs 2 136.0525 <.0001
Fixed effects and test statistics differ between models. A look at marginal
ANOVA tables suggest inference might differ somewhat between models:
anova.lme(m1,type="m")
# numDF denDF F-value p-value
#(Intercept) 1 1789 63384.80 <.0001
#group 445 1.29 0.2893
#x 1 1789 0.05 0.8226
#I(x^2)1 1789 4.02 0.0451
#group:x 4 1789 2.61 0.0341
#group:I(x^2) 4 1789 4.37 0.0016
anova.lme(m2,type="m")
# numDF denDF F-value p-value
#(Intercept) 1 1789 59395.79 <.0001
#group 445 1.33 0.2725
#x 1 1789 0.04 0.8379
#I(x^2)1 1789 2.28 0.1312
#group:x 4 1789 2.09 0.0802
#group:I(x^2) 4 1789 2.81 0.0244
Now, this is all well. But: my colleagues have been running the same data
set using PROC MIXED in SAS and come up with substantially different results
when comparing SAS default covariance structure (variance components) and
AR1. Specifically, there is virtually no change in either test statistics or
fitted values when using AR1 instead of Variance Components in SAS, which
fits the observation that AIC values (in SAS) indicate both covariance
structures fit data equally well.
This is not very satisfactory to me, and I would be interesting to know what
is happening here. Realizing
this might not be the correct forum for this question, I would like to ask
you all if anyone would have any
input as to what is going on here, e.g. am I setting up my model
erroneously, etc.?
N.b. I have no desire to replicate SAS results, but I would most certainly
be interested to know what could possibly explain such a large discrepancy
between the two platforms. Any suggestions greatly welcomed.
(Data are located at:
https://www.dropbox.com/s/hf455aev3teb5e0/data.csv?dl=0)
With all best wishes,
Andreas
--
View this message in context:
http://r.789695.n4.nabble.com/AR1-covariance-structure-for-lme-object-and-R-SAS-differences-in-model-output-tp4703103.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Autoregressive covariance structure for lme object and R/SAS differences in model output
Thanks for your comments. I hard disk meltdown has slowed me down somewhat, but I am now back online. I have reposted to the ME-list, and requested SAS output from my colleagues. I will update the thread again as soon I know more. Thanks, Andreas -- View this message in context: http://r.789695.n4.nabble.com/Autoregressive-covariance-structure-for-lme-object-and-R-SAS-differences-in-model-output-tp4703103p4703396.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] random slope models with lme --> failured to converge
Dear all,
I am working on a data set in which I have sequentially measured egg
temperatures ("eggtemp") in birds incubating in different ambient
temperatures ("treat", sample data set below), "id" is not replicated within
treatment.
id treat eggtemp
1 79 3 30.90166
2 42 3 34.94044
3 10 3 32.69945
4 206 3 36.64127
5 23 3 31.80055
65 3 29.98338
7 24 3 35.72992
8 45 3 30.49584
9 29 3 33.64958
10 78 3 31.37673
11 44 3 32.85873
12 368 3 34.44875
13 79 4 31.24100
14 42 4 34.11634
15 10 4 34.73407
16 206 4 36.20914
17 23 4 34.98061
18 5 4 34.17590
19 24 4 37.71468
20 45 4 35.34765
21 29 4 35.48892
22 78 4 33.26593
23 44 4 34.86981
24 368 4 34.44875
25 79 2 27.33241
26 42 2 30.73269
27 10 2 29.54986
28 206 2 31.78947
29 23 2 29.69114
30 24 2 36.48199
31 45 2 29.76454
32 29 2 30.56510
33 78 2 27.71468
For this data, I want to construct a model with a random intercept and slope
to compare with a model containing only the random intercept to check
whether or not different individuals respond differently to treatment.
> temp2.lme<-lme(eggtemp~treat,random=~treat|id,data=temp3.df)
However, I get the below error, which I suspect might be because individuals
are not replicated within treatments.
> nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (9)
If I specify treatment as a factor (which it really is not, just different
ambient temperatures), the model seems to converge. Similarly, the model
converges if I exclude one of the treatment levels. However, this does not
feel very satisfactorily.
So, how do I move on from here? Any tips and hints are much appreciated. I
have tried to include the random slope only, which also works, but the
inference for such a model would be quite different and not really what I am
after.
> temp3.lme<-lme(eggtemp~treat,random=~treat-1|id,data=temp3.df)
Kind regards,
Andreas Nord
Lund University
Sweden
--
View this message in context:
http://n4.nabble.com/random-slope-models-with-lme-failured-to-converge-tp1469575p1469575.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Problems with lme random slope+intercept model
Dear all, I'm trying to fit a model on ecological data in which I have measured a few biotic and abiotic factors over the course of a few days in several individuals. Specifically, I'm interested in modelling y ~ x1, with x2, x3, and 'factor' as independent variables. Because data suggests both slope and intercept (for y ~x1) might differ between individuals, I'd want to compare model fit for a saturated model with random intercept only, against that of a model with random slope + intercept. Data are available in full from this link: https://www.dropbox.com/s/mzk8utvgkzp4rtr/data.txt The random intercept model seems to function appropriately: data<-subset(data,data$id!='id225' & data$id!='id237' & data$id!='id233') m1.lme<-with(data,lme(y~x1+x2+x3+factor,random=~1|id,na.action=na.omit)) However, fitting the random slope+intercept model produces an error message I can't quite make sense of. m2.lme<-with(data,lme(y~x1+x2+x3+factor,random=~1+y|id,na.action=na.omit)) #Error in chol.default((value + t(value))/2) : # the leading minor of order 2 is not positive definite I also tried fitting the same model with a diagonal covariance structure, which resulted in convergence failure. m3.lme<-with(data,lme(y~x1+x2+x3+factor,random=reStruct(object=~1+y|id,pdClass="pdDiag"),na.action=na.omit)) #Error in lme.formula(y ~ x1 + x2 + x3, random = reStruct(object = ~y | : #nlminb problem, convergence error code = 1 #message = false convergence (8) However, changing lmeControl gets this model to run, but I can't make sense of the estimates for fixed effects, suggesting the model might be biased. In addition, I'm not sure how changing lmeControl changes model interpretation. Perhaps someone could fill me in on this? m4.lme<-with(data,lme(y~x1+x2+x3+factor,random=reStruct(object=~y|id,pdClass="pdDiag"),na.action=na.omit, control=lmeControl(opt = "optim"))) Any hints on how to proceed from this would be greatly appreciated. Best, and thanks, Andreas -- View this message in context: http://r.789695.n4.nabble.com/Problems-with-lme-random-slope-intercept-model-tp4679104.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problems with lme random slope+intercept model
Dear Ben, Thank you, I agree. I have forwarded this to the r-sig-mixed-models list. Best, Andreas -- View this message in context: http://r.789695.n4.nabble.com/Problems-with-lme-random-slope-intercept-model-tp4679104p4679168.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Crossed random effects in lme
Dear all, I am trying to fit a model with crossed random effects using lme. In this experiment, I have been measuring oxygen consumption (mlmin) in bird nestlings, originating from three different treatments (treat), in a respirometer with 7 different channels (ch). I have also measured body mass (mass) for these birds. id nesttreat yearmlmin massch hack 1EP5171117 H 20081.401719138 10.74 2008:17 1EP5170917 H 20081.257163112 9.7 5 2008:17 1EP5171617 H 20081.050170714 10.26 2008:17 1EP5171217 H 20081.330495314 9.6 7 2008:17 1EP51791687 M 20081.07625708 9.7 3 2008:687 1EP51772823 H 20081.336820232 10.24 2008:823 1EP51778613 L 20081.300814516 10.75 2008:613 1EP52336207 M 20081.071775936 10.73 2008:207 1EP52403808 H 20081.142389688 10.35 2008:808 1ER17603838 M 20090.984225217 9.6 3 2009:838 1ER17607838 M 20091.045058894 9.3 4 2009:838 1ER17600247 L 20091.047603048 9.2 5 2009:247 1ER17299247 L 20090.974569658 9.2 6 2009:247 1ER17292617 H 20091.271260094 10.57 2009:617 1ER172067009M 20091.074791644 10.72 2009:7009 1ER17221730 H 20091.423266177 10.24 2009:730 1ER17275863 L 20091.433076022 10.74 2009:863 1ER17277863 L 20091.165236024 9.7 5 2009:863 1ER17283863 L 20091.139311895 10.46 2009:863 1ER17280863 L 20091.056161196 10.47 2009:863 CK59991 690 H 20100.994878996 9.5 2 2010:690 CK59806 161 M 20101.070052025 9.7 6 2010:161 CK59859 545 M 20101.456680579 9.9 4 2010:545 CK59862 545 M 20101.350698793 9.9 5 2010:545 CK59871 223 L 20100.830582186 8.3 6 2010:223 CK59868 223 L 20100.776241825 8 7 2010:223 CL77343 365 M 20101.352454484 10.34 2010:365 CL77338 365 M 20101.327691628 9.6 5 2010:365 CL77356 191 H 20101.212796979 11.31 2010:191 CL77361 191 H 20100.882307732 11.42 2010:191 CL77355 191 H 20101.137097586 10.93 2010:191 I want to include both nesting attempt (hack) and respirometer channel (ch) as random factors in a model trying to explain variation in oxygen consumption. From Pinheiro & Bates (2000), I've gathered that this model could be fit making use of pdBlocked and pdIdent, so I've tried fitting the below model: m1.bmr<-with(bmred.df,lme(mlmin~treat*year+massout,random=pdBlocked(list(pdIdent(~hack-1),pdIdent(~ch-1))) )) However, my model fails with the following error message: Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups I would much appreciate any input on this! Kind regards, Andreas Nord Sweden -- View this message in context: http://r.789695.n4.nabble.com/Crossed-random-effects-in-lme-tp3000101p3000101.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Problems with glht function for lme object
Dear all,
I'm trying to make multiple comparisons for an lme-object. The data is for
an experiment on parental work load in birds, in which adults at different
sites were induced to work at one of three levels ('treat'; H, M, L). The
response is 'feedings', which is a quantitative measure of nest provisioning
per parent per chick per hour. Site is included as a random effect (one
male/female pair per site).
My final model takes the following form:
feedings ~ treat + year + data^2, random = ~1|site,data=feed.df
For this model, I would like to do multiple comparisons on 'treat', using
the multcomp package:
summary(glht(m4.feed,linfct=mcp(treat="Tukey")))
However, this does not work, and I get the below error message.
Error in if (is.null(pkg) | pkg == "nlme") terms(formula(x)) else slot(x, :
argument is of length zero
Error in factor_contrasts(model) :
no ‘model.matrix’ method for ‘model’ found!
I suspect this might have quite a straightforward solution, but I'm stuck at
this point.
Any help would be most appreciated. Sample data below.
Kind regards,
Andreas Nord
Sweden
==
feedings sex site treat year date^2
1.8877888 M 838 H 2009 81
1.9102787 M 247 H 2009 81
1.4647229 M 674 H 2010 121
1.4160590 M 7009 M 2009 144
1.3106749 M 863 M 2010 196
1.2718121 M 61 M 2009 225
1.2799263 M 729 L 2009 256
1.5829564 M 629 L 2009 256
1.4847251 M 299 L 2010 324
1.2463151 M 569 L 2010 324
2.1694169 F 838 H 2009 81
1.5966899 F 247 H 2009 81
2.4136983 F 674 H 2010 121
1.7784873 F 7009 M 2009 144
1.6681317 F 863 M 2010 196
2.3691275 F 61 M 2009 225
2.0672192 F 729 L 2009 256
1.6389902 F 629 L 2009 256
0.9307536 F 299 L 2010 324
1.6786767 F 569 L 2010 324
==
--
View this message in context:
http://r.789695.n4.nabble.com/Problems-with-glht-function-for-lme-object-tp3179128p3179128.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problems with glht function for lme object
Dear all, Thanks for your input. It works fine now. All it took was for me to tidy up the workspace. Simple solution, that I should have considered earlier. Regards, Andreas -- View this message in context: http://r.789695.n4.nabble.com/Problems-with-glht-function-for-lme-object-tp3179128p3206686.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] glht for lme object with significant interaction term
Dear all, I'm working on some data from an experiment on the breeding behavior of birds. In short, I have been measuring how the time spent on performing a certain task (variable 'mean_on_active') differs over time (variable 'day', 2 levels) across three experimental categories (variable 'treat'; levels 'C', 'R', 'E'). The model shows a significant interaction between treatment and day. To have a closer look at this, I would like to do multiple comparisons for treatment levels within day (i.e. across treatments for each day in turn). I have browsed the forum for a way of doing this using the glht function, but haven't found a good solution to my problem. Any hints on how to proceed, or on other methods that might be (more) appropriate? Final model output below, data attached. http://r.789695.n4.nabble.com/file/n4097865/sub.data.txt sub.data.txt Kind regards, Andreas ##Final model m.final<-lme(mean.on.active~treat+day+treat:day,random=1|id,na.action=na.omit) summary(m.final) Linear mixed-effects model fit by REML Data: NULL AIC BIClogLik 1051.779 1075.757 -517.8896 Random effects: Formula: ~1 | id (Intercept) Residual StdDev:6.276309 5.473167 Fixed effects: mean_on_active ~ treat * day Value Std.Error DF t-value p-value (Intercept) 24.987123 1.792599 75 13.939049 0. treatE15.661119 2.327329 73 6.729225 0. treatR 1.133678 2.228713 73 0.508670 0.6125 day6-7-1.068160 1.758899 73 -0.607289 0.5455 treatE:day6-7 -6.650690 2.335567 73 -2.847570 0.0057 treatR:day6-7 -1.495927 2.273071 73 -0.658108 0.5125 Correlation: (Intr) treatE treatR day6-7 tE:6-7 treatE-0.736 treatR-0.747 0.572 day6-7-0.482 0.372 0.389 treatE:day6-7 0.361 -0.522 -0.292 -0.753 treatR:day6-7 0.370 -0.281 -0.524 -0.774 0.582 Standardized Within-Group Residuals: Min Q1Med Q3Max -2.5774244 -0.4585991 -0.1360731 0.3902143 3.9752257 Number of Observations: 154 Number of Groups: 76 -- View this message in context: http://r.789695.n4.nabble.com/glht-for-lme-object-with-significant-interaction-term-tp4097865p4097865.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
