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

Sebastián Daza sebastian.daza at gmail.com
Tue Apr 12 22:47:01 CEST 2011 

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

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