[R-sig-ME] using weights=varIdent in lme: how to change reference level of grouping factor?
Karl Ove Hufthammer
karl at huftis.org
Wed Jul 17 19:12:58 CEST 2013
Ben Bolker wrote:
>>I am fitting a model using the lme function. I am using
>> weights=varIdent to fit a variance model with different variances
>> for each level of a stratification variable.Can someone please tell
>> me how to chance the reference level of the stratification variable?
>> I have tried changing the reference level of the factor, and tried
>> re-ordering the data frame so that the first observation has the
>> level that I want to use as the reference level, but neither have
>> worked. Emma
>
> Good question. I would have guessed that re-ordering the factor
> would do the trick, but I guess not.
Is the following problem related? When using *initalized* variance
functions, strange things happen:
$ library(nlme)
$ Orth2 <- transform(Orthodont,
+ Sex=factor(Sex,levels=c("Female","Male")))
$ Orth3=Orthodont[nrow(Orthodont):1,]
$ var1=Initialize(varIdent(form=~1|Sex), data=Orthodont)
$ var2=Initialize(varIdent(form=~1|Sex), data=Orth2)
$ var3=Initialize(varIdent(form=~1|Sex), data=Orth3)
$ var1
Variance function structure of class varIdent representing
Male Female
1 1
$ var2
Variance function structure of class varIdent representing
Male Female
1 1
$ var3
Variance function structure of class varIdent representing
Female Male
1 1
So reordering the levels doesn’t change the reference level (var2), but
reordering the data seems to (var3). That’s nice. But when we actually try
to fit the models, we get nonsensical results:
$ fm1=lme(distance ~ age, random =~age|Subject,
+ weights=var1, data=Orthodont)
$ fm2=lme(distance ~ age, random =~age|Subject, weights=var2, data=Orth2)
$ fm3=lme(distance ~ age, random =~age|Subject, weights=var3, data=Orth3)
$ fm1$modelStruct$varStruct
Variance function structure of class varIdent representing
Male Female
1.000000 0.404098
$ fm2$modelStruct$varStruct
Variance function structure of class varIdent representing
Male Female
1.000000 0.404098
$ fm3$modelStruct$varStruct
Variance function structure of class varIdent representing
Female Male
1.0000000 0.4315485
This result seems completely wrong. Shouldn’t the result for Male be ~2.4?
The estimates are identical for the first two models, but *slightly*
different for the third one:
$ fixef(fm1)
(Intercept) age
17.2983171 0.6068598
$ fixef(fm2)
(Intercept) age
17.2983171 0.6068598
$ fixef(fm3)
(Intercept) age
17.3879588 0.6046554
But the most surprising thing (for me) is that the model fitting doesn’t
seem to care about the value for the stratifying variable at all. For
example, running
Orth4=transform(Orthodont, Sex="test")
And fitting the model gives the same results as for fm3. I thought maybe the
reason was the lme used the values from the variance function instead of
from the data frame, but even if I completely reorder the data frame,
breaking any implicit match, I get the same results as for fm3:
$ Orth4=transform(Orthodont, Sex="test")
$ Orth4=Orth4[sample(nrow(Orth4)),]
$ fm4=lme(distance ~ age, random =~age|Subject, weights=var3,
+ data=Orth4)
$ fm4$modelStruct$varStruct
Variance function structure of class varIdent representing
Female Male
1.0000000 0.4315485
$ fixef(fm4)
(Intercept) age
17.3879588 0.6046554
Perhaps the reason could be that the differences in variances between the
strata are so small that they don’t have any effect on estimates? However,
fitting a model with *equal* within-group variances does give somewhat
different estimates:
$ fm5=lme(distance ~ age, random =~age|Subject, data=Orth4)
$ fixef(fm5)
(Intercept) age
16.7611111 0.6601852
Version info:
$ sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-suse-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=nn_NO.UTF-8 LC_NUMERIC=C
LC_TIME=nn_NO.UTF-8 LC_COLLATE=nn_NO.UTF-8
LC_MONETARY=nn_NO.UTF-8
[6] LC_MESSAGES=nn_NO.UTF-8 LC_PAPER=C
LC_NAME=C LC_ADDRESS=C
LC_TELEPHONE=C
[11] LC_MEASUREMENT=nn_NO.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] nlme_3.1-110
loaded via a namespace (and not attached):
[1] grid_3.0.1 lattice_0.20-15 tools_3.0.1
--
Karl Ove Hufthammer
E-mail/Jabber: karl at huftis.org
More information about the R-sig-mixed-models
mailing list