[R-sig-ME] level 1 variance-covariance structure

Tue Apr 12 17:21:24 CEST 2011

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
> 