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

Tue May 19 11:05:48 CEST 2015

Dear Asher,

I would focus on the mathematical expressions used for both the R and SPSS
models. Then it should be straightforward to see if both models fit the
same thing. And if so which parameters have the same interpretation.

You can find the expressions for nlme in Pinheiro and Bates (2000). I can't
help you with the SPSS part.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2015-05-19 10:55 GMT+02:00 Asher Strauss <asher.strauss op gmail.com>:

> 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_Beyond_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
>
> On Mon, May 18, 2015 at 9:01 PM, Daniel Wright <Daniel.Wright op act.org>
> wrote:
>
> > It may also be convenient to use the gls function in nlme.
> >
> > This is used in
> >
> https://www.researchgate.net/publication/23134911_Multilevel_modelling_Beyond_the_basic_applications
> > and in http://www.ats.ucla.edu/stat/r/examples/alda/ch7.htm which is the
> > UCLA page for Singer and Willet's wonderful book.
> >
> >
> > -----Original Message-----
> > From: R-sig-mixed-models [mailto:
> r-sig-mixed-models-bounces op r-project.org]
> > On Behalf Of Asher Strauss
> > Sent: Saturday, May 16, 2015 1:19 AM
> > To: r-sig-mixed-models op r-project.org
> > Subject: [R-sig-ME] estimating AR1 parameters of level one error using
> lme
> >
> > Hi all,
> > I am new to using R for mixed effects (have been using SPSS until now),
> > please for give me if this is a trivial question.
> >
> > I am trying to understand how to estimate AR1 parameters of level one
> > error using lme (I have understood that specifying level one
> > variance-covariance matrix is not easily possible in lmer, is this
> true?).
> >
> > When using SPSS one estimates two parameters: AR1 diagonal and AR1 rho. I
> > am searching for an equivalent command in R.
> >
> > 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 out put of:
> >
> > Correlation Structure: ARMA(1,0)
> >  Formula: ~1 + Time | Subject
> >  Parameter estimate(s):
> >      Phi1
> > 0.4334469
> >
> > is Phi1 equivalent to Rho?  I do not believe so, since when estimating
> AR1
> > diagonal and AR1 rho using SPSS I received 1.349 and -0.942 respectively.
> >
> > here is the SPSS syntax I am using:
> > 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).
> >
> > Thank you!
> > Asher
> >
> >         [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models op r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]