Re: [R] AR1 covariance structure for lme object and R/SAS differences in model output
anord biol.lu.se> writes: > [snip snip] > 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) Could you repost this on r-sig-mixed-mod...@r-project.org ? It would be useful to see some of the comparisons (fixed effects, RE variance-covariance, correlation parameter) between SAS and R; are they actually fitting the same model? (e.g., does SAS allow for covariance between the slope and intercept random effects?) Maybe they're getting the same estimates but computing df/p-values in different ways? I thought I would try this with orthogonal polynomials in case some of the fits were unstable ... data.df <- read.csv2("ar1_data.csv") library("nlme") m1 <- lme(y~group*x+group*I(x^2),random=~x|id, data=data.df,na.action=na.omit,method="ML") ## use method="ML" so we can compare orthog. and non-orthog. polynomials m1B <- update(m1,.~group*poly(x,2)) m2 <- update(m1,corr=corCAR1(form=~x|id)) m2B <- update(m1B,corr=corCAR1(form=~x|id)) AIC(m1,m1B,m2,m2B) __ 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] AR1 covariance structure for lme object and R/SAS differences in model output
> On 11 Feb 2015, at 16:57 , anord wrote: > > 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) > Hmm, does SAS include a nugget effect perchance? At any rate, showing the SAS output (or some of it) might make it easier for someone to help. Also, R-sig-ME is a better choice for discussions of mixed effects models. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.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] 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.