level 1 variance-covariance structure
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