[R-sig-ME] level 1 variance-covariance structure
Sebastián Daza
sebastian.daza at gmail.com
Thu Apr 14 00:28:46 CEST 2011
Finally, I had to run the models using HLM, where I could get results
(without iteration errors) for a model like this:
attit ~ age13 , data, random= ~ age13 | id, correlation = corAR1(form =
~ wave | id))
I don't know why I could run this kind of models in HLM and I couldn't
do it using R (lme). It would be good to know if it is a limitation of
the R package or an over-parametrization of the model... or it is
related to a specific characteristic of the dataset. I don't have any clue.
Thank you!
On 4/13/2011 3:29 AM, ONKELINX, Thierry wrote:
> 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
>>
--
Sebastián Daza
sebastian.daza at gmail.com
More information about the R-sig-mixed-models
mailing list