[R-sig-ME] estimating AR1 parameters of level one error using lme

Viechtbauer Wolfgang (STAT) wolfgang.viechtbauer at maastrichtuniversity.nl
Tue May 19 14:55:28 CEST 2015


First of all, double, no, tripple check that the data are exactly the same in R and in SPSS. I have seen too many cases of "Why are my results in R different than in <insert other software package here>?" that could be traced down to export/import issues/mistakes (e.g., things like -999 that are declared to be missing values in SPSS, but that are read in as is into R).

Assuming that the data are in fact the same -- yes, what is labeled 'Phi1' in the output from lme() is indeed the estimate of the autocorrelation parameter. However, there is already one thing that is different between your two models. In R, you used 'random=~1+Time|Subject' (you could also just write 'random=~Time|Subject' -- that does the same thing). That will add random intercepts and slopes to the model that are allowed to be correlated. In SPSS, you used '/RANDOM=INTERCEPT Time | SUBJECT(Subject)'. If I am not mistaken, the default is COVTYPE(VC), which adds random intercepts and slopes to the model, which are assumed to be independent. So, your models are not identical to begin with.

I would suggest to first fit the models without the AR(1) structure and make sure you get essentially the same results. Then you can add the AR(1) part and cross-check the results.

Best,
Wolfgang

--    
Wolfgang Viechtbauer, Ph.D., Statistician    
Department of Psychiatry and Neuropsychology    
School for Mental Health and Neuroscience    
Faculty of Health, Medicine, and Life Sciences    
Maastricht University, P.O. Box 616 (VIJV1)    
6200 MD Maastricht, The Netherlands    
+31 (43) 388-4170 | http://www.wvbauer.com    

> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-
> project.org] On Behalf Of Asher Strauss
> Sent: Tuesday, May 19, 2015 10:56
> To: Daniel Wright
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] estimating AR1 parameters of level one error
> using lme
> 
> Hi Daniel
> and thank you very much for your reply!
> 
> I went over the documents you sent me, though couldn't find the sdtalt
> package :(
> 
> I think using  lme is more suitable for my purposes since it is important
> for me include random effects in my model.
> 
> I am still stuck with understanding the right code in R to specify a
> homogeneous variance structure with correlations structured as AR1 (model
> 13 in your paper
> https://www.researchgate.net/publication/23134911_Multilevel_modelling_Be
> yond_the_basic_applications
> ).
> As I wrote in my previous message specifying this model in R (if I got
> the
> code right - I assume not...) and in SPSS produces completely different
> results.
> 
> I will be grateful if you or anybody else can point me towards a
> solution.
> 
> I am repeating here my previous example to make my difficulty easier to
> address:
> 
> for example using the Glucose data from the nlme package:
> 
> data(Glucose)
> 
> fGlucose<-filter(Glucose,Meal=="10am")
> 
> summary(
> lme(
> fixed=conc~Time,
> random=~1+Time|Subject,
> method="REML",
> data=fGlucose,
> na.action="na.omit",
> correlation=corAR1(form=~1+Time|Subject))
> )
> 
> 
> I get an output of:
> 
> Correlation Structure: ARMA(1,0)
>  Formula: ~1 + Time | Subject
>  Parameter estimate(s):
>      Phi1
> 0.4334469
> 
> I am assuming  Phi1 is equivalent to Rho (am I right?). What makes me
> suspect I am wrong, and confuses me, is that when estimating AR1
> diagonal (homogeneous
> level one error variance) and AR1 rho (the auto regression parameter)
>  using SPSS I received 1.349 and -0.942 respectively.
> 
> here is the SPSS syntax I am using  on the same data set (Glucose):
> 
> COMPUTE filter_$=(Meal="10am").
> VARIABLE LABELS filter_$ 'Meal="10am" (FILTER)'.
> VALUE LABELS filter_$ 0 'Not Selected' 1 'Selected'.
> FORMATS filter_$ (f1.0).
> FILTER BY filter_$.
> EXECUTE.
> 
> MIXED conc WITH Time
>  /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
> SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE)
> PCONVERGE(0.000001, ABSOLUTE)
>   /FIXED=INTERCEPT Time | SSTYPE(3)
>   /METHOD=REML
>   /PRINT=SOLUTION TESTCOV
>   /RANDOM=INTERCEPT Time | SUBJECT(Subject)
>   /REPEATED=Time | SUBJECT(Subject)COVTYPE(AR1).
> 
> It would be really great if you or anybody else can assist me here!
> Thank you!
> 
> Best
> Asher Strauss



More information about the R-sig-mixed-models mailing list