[R-sig-ME] nlme package: changing reference for varIdent parameter estimates in summary.gls
Kornak, John
john.kornak at ucsf.edu
Wed Dec 24 23:34:09 CET 2014
Dear R mixed-effects modeling experts,
I am running gls models with heterogeneous group variances using the glm function in the nlme package with varIdent weights. I am interested in controlling the baseline level for the group variance function standard deviation estimates by applying the relevel function to the group variable. When I use relevel, the baseline level in the coefficient table changes as expected, but the baseline level does not change for the variance function table; for example, see the fitted models gls1 vs. gls2 in the contrived example below. Does anyone have a suggested hack as to how I can change the baseline level for the variance function output please? One hack I have made work is to change the order of the rows in the data (see final part of example), but that seems rather ugly and will get painful for the real data I am working with for which I need to reset many different reference levels.
I have google searched / checked help pages for other solutions without success.
I am running R version 3.1.0 on an iMac OSX v. 10.9.5
Thank you in advance and happy holidays!
John Kornak
> library(nlme)
> group <- factor(c(rep("no",20),rep("yes",20)))
> set.seed(2)
> outcome <- c(rnorm(20,0,2),rnorm(20,5,4))
> dataTest <- data.frame(outcome,group)
# Original model fit before releveling
> gls1 <- gls(outcome ~ group, weights=varIdent(form = ~1|group), data=dataTest)
> summary(gls1)
— snip —
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | group
Parameter estimates:
no yes
1.00000 2.23034
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.390922 0.4734001 0.825775 0.4141
groupyes 4.607951 1.1571140 3.982279 0.0003
—snip —
Residual standard error: 2.11711
Degrees of freedom: 40 total; 38 residual
# relevel the group so that “yes” is the reference
> dataTest$group <- relevel(dataTest$group,"yes")
> gls2 <- gls(outcome ~ group, weights=varIdent(form = ~1|group), data=dataTest)
> summary(gls2)
— snip —
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | group
Parameter estimates:
no yes
1.00000 2.23034 ### “no" is still the reference group here for the variance function
Coefficients:
Value Std.Error t-value p-value
(Intercept) 4.998873 1.055843 4.734484 0e+00
groupno -4.607951 1.157114 -3.982279 3e-04 # “yes” has become the reference for the coefficients
— snip —
Residual standard error: 2.11711
Degrees of freedom: 40 total; 38 residual
>
# Reorganize the data so that the group == “yes" entries come first
> dataTest2 <- dataTest[c(21:40,1:20), ]
> gls3 <- gls(outcome ~ group, weights=varIdent(value = c(yes = 0.5), form = ~1|group), data=dataTest2)
> summary(gls3)
— snip —
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | group
Parameter estimates:
yes no
1.0000000 0.4483626 ## now I have the desired “yes” as reference level
Coefficients:
Value Std.Error t-value p-value
(Intercept) 4.998873 1.055843 4.734487 0e+00
groupno -4.607951 1.157113 -3.982281 3e-04 # and “yes” also reference level here
— snip —
Residual standard error: 4.721872
Degrees of freedom: 40 total; 38 residual
>
---------------------------------------------------
John Kornak, PhD
Associate Professor in Residence
Department of Epidemiology and Biostatistics
University of California, San Francisco
Mission Hall: Global Health & Clinical Sciences Building
550 16th St, 2nd floor, Box #0560
San Francisco, CA 94158-2549
Tel: 415-514-8028
Fax: 415-514-8150
Email: john.kornak at ucsf.edu<mailto:john.kornak at ucsf.edu>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list