[R-sig-ME] level 1 variance-covariance structure
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Wed Apr 13 10:29:38 CEST 2011
There is no auto-correlation left AFTER the fixed and random effects are taken into account. So you probably will have to choose between the models below.
m3a <- lme(attit ~ age13 , data, random= ~ 0 + factor(age13)| id)
m3b <- lme(attit ~ age13 , data, random= ~ 1| id, correlation = corAR1(form = ~ age13 | id))
----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be
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
> -----Oorspronkelijk bericht-----
> Van: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens
> Sebastián Daza
> Verzonden: dinsdag 12 april 2011 22:47
> Aan: R-SIG-Mixed-Models at r-project.org
> Onderwerp: Re: [R-sig-ME] level 1 variance-covariance structure
>
> Thierry,
> I can run this model... but what does it mean?
> The correlation structure that I get is:
>
> Correlation Structure: ARMA(1,0)
> Formula: ~age13 | id
> Parameter estimate(s):
> Phi1
> 0
>
> What does zero mean? I would expect get some positive number there...
>
> m3a <- lme(attit ~ 1 + age13 , data, random= ~ 0 +
> factor(age13)| id, correlation = corAR1(form = ~ age13 | id))
> > summary(m3a)
> Linear mixed-effects model fit by REML
> Data: data
> AIC BIC logLik
> -324.2096 -229.5528 181.1048
>
> Random effects:
> Formula: ~0 + factor(age13) | id
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> factor(age13)-2 0.17219431 f(13)-2 f(13)-1 f(13)0 f(13)1
> factor(age13)-1 0.19789254 0.493
> factor(age13)0 0.25942941 0.425 0.544
> factor(age13)1 0.28354459 0.442 0.442 0.723
> factor(age13)2 0.29097081 0.498 0.474 0.639 0.808
> Residual 0.07457025
>
> Correlation Structure: ARMA(1,0)
> Formula: ~age13 | id
> Parameter estimate(s):
> Phi1
> 0
> Fixed effects: attit ~ 1 + age13
> Value Std.Error DF t-value p-value
> (Intercept) 0.3210558 0.012832840 839 25.01829 0
> age13 0.0593529 0.004716984 839 12.58282 0
> Correlation:
> (Intr)
> age13 0.504
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -1.46371874 -0.27170442 -0.04080688 0.26239551 1.69883907
>
> Number of Observations: 1079
> Number of Groups: 239
>
> On 4/12/2011 10:21 AM, ONKELINX, Thierry wrote:
> > Dear Sebastian,
> >
> > The model below works fine on my computer.
> >
> > m3a<- lme(attit ~ 1 + age13 , data=dataset, random= ~
> > 0+factor(age13)| id, correlation = corAR1(form = ~ age13 | id))
> >
> >
> > Best regards,
> >
> > Thierry
> >
> ----------------------------------------------------------------------
> > ------
> > ir. Thierry Onkelinx
> > Instituut voor natuur- en bosonderzoek team Biometrie&
> Kwaliteitszorg
> > Gaverstraat 4 9500 Geraardsbergen Belgium
> >
> > Research Institute for Nature and Forest team Biometrics& Quality
> > Assurance Gaverstraat 4 9500 Geraardsbergen Belgium
> >
> > tel. + 32 54/436 185
> > Thierry.Onkelinx at inbo.be
> > www.inbo.be
> >
> > 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
> >
> >
> >> -----Oorspronkelijk bericht-----
> >> Van: Sebastián Daza [mailto:sebastian.daza at gmail.com]
> >> Verzonden: dinsdag 12 april 2011 15:43
> >> Aan: ONKELINX, Thierry
> >> CC: R-SIG-Mixed-Models at r-project.org
> >> Onderwerp: Re: [R-sig-ME] level 1 variance-covariance structure
> >>
> >> Thank you for your reply Thierry...
> >> Increasing the number of iterations doesn't work:
> >>
> >> m3a<- lme(attit ~ 1 + age13 , data=data, random= ~ age13 | id,
> >> correlation = corAR1(, form = ~ ind | id),
> >> control=list(maxIter=1000, msMaxIter=1000, niterEM=1000))
> >>
> >> Error in lme.formula(attit ~ 1 + age13, data = data, random =
> >> ~age13 | :
> >> nlminb problem, convergence error code = 1
> >> message = function evaluation limit reached without
> convergence
> >> (9)
> >>
> >> I have attached my database. I don't know if it is a problem of my
> >> model or a limitation of lme function.
> >>
> >> The best!
> >> Sebastian.
> >>
> >> On 4/12/2011 6:25 AM, ONKELINX, Thierry wrote:
> >>> Dear Sebastian,
> >>>
> >>> You don't need to create dummy variables your selve.
> >>>
> >>> You can write m2a<- lme(attit ~ 1 + age13 , data=data,
> >> random= ~ 0 + ind1+ ind2+ ind3+ ind4+ ind5 | id, method="REML") as
> >>>
> >>> m2a<- lme(attit ~ 1 + age13 , data=data, random= ~ 0 +
> >> factor(ind) | id, method="REML")
> >>>
> >>> Or if ind is an indicator for age13:
> >>>
> >>> m2a<- lme(attit ~ 1 + age13 , data=data, random= ~ 0 +
> >> factor(age13) | id, method="REML")
> >>>
> >>> Have a look at lmeControl() to increase the number of iterations.
> >>>
> >>> Best regards,
> >>>
> >>> Thierry
> >>>
> >>>
> >> --------------------------------------------------------------
> >> --------------
> >>> ir. Thierry Onkelinx
> >>> Instituut voor natuur- en bosonderzoek
> >>> team Biometrie& Kwaliteitszorg
> >>> Gaverstraat 4
> >>> 9500 Geraardsbergen
> >>> Belgium
> >>>
> >>> Research Institute for Nature and Forest
> >>> team Biometrics& Quality Assurance
> >>> Gaverstraat 4
> >>> 9500 Geraardsbergen
> >>> Belgium
> >>>
> >>> tel. + 32 54/436 185
> >>> Thierry.Onkelinx at inbo.be
> >>> www.inbo.be
> >>>
> >>> 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
> >>>
> >>>
> >>>> -----Oorspronkelijk bericht-----
> >>>> Van: r-sig-mixed-models-bounces at r-project.org
> >>>> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens
> Sebastián
> >>>> Daza
> >>>> Verzonden: maandag 11 april 2011 18:44
> >>>> Aan: R-SIG-Mixed-Models at r-project.org
> >>>> Onderwerp: [R-sig-ME] level 1 variance-covariance structure
> >>>>
> >>>> Hi everyone,
> >>>> I am trying to reproduce some results models from HLM (HMLM) to
> >>>> contrast different specifications of level 1
> variance-covariance,
> >>>> but I get convergence errors. I would like to know if
> there are any
> >>>> problems with my model
> >> specification...
> >>>>
> >>>>
> >>>> # database structure
> >>>>
> >>>> head(data[,c(1,2,6, 9:13,17)])
> >>>> id attit age13 ind1 ind2 ind3 ind4 ind5 ind
> >>>> 1 3 0.11 -2 1 0 0 0 0 1
> >>>> 2 3 0.20 -1 0 1 0 0 0 2
> >>>> 3 3 0.00 0 0 0 1 0 0 3
> >>>> 4 3 0.00 1 0 0 0 1 0 4
> >>>> 5 3 0.11 2 0 0 0 0 1 5
> >>>> 6 8 0.29 -2 1 0 0 0 0 1
> >>>>
> >>>> # attit is a deviant measure and ind variables indicate
> different
> >>>> waves # following some examples of snijders and bosker's book, I
> >>>> get the unrestricted model:
> >>>>
> >>>> > m2a<- lme(attit ~ 1 + age13 , data=data, random= ~ 0 +
> >>>> ind1+ ind2+
> >>>> ind3+ ind4+ ind5 | id, method="REML")
> >>>>
> >>>> > summary(m2a)
> >>>> Linear mixed-effects model fit by REML
> >>>> Data: data
> >>>> AIC BIC logLik
> >>>> -326.2096 -236.5348 181.1048
> >>>>
> >>>> Random effects:
> >>>> Formula: ~0 + ind1 + ind2 + ind3 + ind4 + ind5 | id
> >>>> Structure: General positive-definite, Log-Cholesky
> >> parametrization
> >>>> StdDev Corr
> >>>> ind1 0.17219431 ind1 ind2 ind3 ind4
> >>>> ind2 0.19789253 0.493
> >>>> ind3 0.25942942 0.425 0.544
> >>>> ind4 0.28354459 0.442 0.442 0.723
> >>>> ind5 0.29097082 0.498 0.474 0.639 0.808
> >>>> Residual 0.07457025
> >>>>
> >>>> Fixed effects: attit ~ 1 + age13
> >>>> Value Std.Error DF t-value p-value
> >>>> (Intercept) 0.3210558 0.012832840 839 25.01829 0
> >>>> age13 0.0593529 0.004716984 839 12.58282 0
> >>>> Correlation:
> >>>> (Intr)
> >>>> age13 0.504
> >>>>
> >>>> Standardized Within-Group Residuals:
> >>>> Min Q1 Med Q3 Max
> >>>> -1.46371871 -0.27170442 -0.04080686 0.26239553 1.69883910
> >>>>
> >>>> Number of Observations: 1079
> >>>> Number of Groups: 239
> >>>>
> >>>> # variance-covariance matrix
> >>>>
> >>>> > extract.lme.cov2(m2a,data)$V[[6]]
> >>>> 25 26 27 28 29
> >>>> 25 0.03521160 0.01681647 0.01899029 0.02159300 0.02494013
> >>>> 26 0.01681647 0.04472218 0.02793174 0.02481343 0.02727012
> >>>> 27 0.01899029 0.02793174 0.07286434 0.05318967 0.04823107
> >>>> 28 0.02159300 0.02481343 0.05318967 0.08595826 0.06667139
> >>>> 29 0.02494013 0.02727012 0.04823107 0.06667139 0.09022474
> >>>>
> >>>> # I get the same results than unrestricted model in HLM
> >>>>
> >>>> # When I try to get the same unrestricted model using "corStruc"
> >>>> commands in lme, I get a convergence problem. Am I
> reproducing the
> >>>> model m2a?
> >>>>
> >>>> > m2b<- lme(attit ~ 1 + age13 , data=data, random= ~ age13
> >>>> | id, correlation = corSymm(, form = ~ ind | id)) Error in
> >>>> lme.formula(attit ~ 1 + age13, data = data, random = ~age13 | :
> >>>> nlminb problem, convergence error code = 1
> >>>> message = iteration limit reached without convergence (9)
> >>>>
> >>>> # When I try to get an autoregressive model, I get again a
> >>>> convergence problem.
> >>>>
> >>>> > m3a<- lme(attit ~ 1 + age13 , data=data, random= ~ age13
> >>>> | id, correlation = corAR1(, form = ~ ind | id)) Error in
> >>>> lme.formula(attit ~ 1 + age13, data = data, random = ~age13 | :
> >>>> nlminb problem, convergence error code = 1
> >>>> message = iteration limit reached without convergence (9)
> >>>>
> >>>> Does anyone know how I can solve this?
> >>>> Thank you in advance.
> >>>>
> >>>> --
> >>>> Sebastián Daza
> >>>> sebastian.daza at gmail.com
> >>>>
> >>>> _______________________________________________
> >>>> R-sig-mixed-models at r-project.org mailing list
> >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>>>
> >>
> >> --
> >> Sebastián Daza
> >> sebastian.daza at gmail.com
> >>
>
> --
> Sebastián Daza
> sebastian.daza at gmail.com
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list