[R-sig-ME] nlme package: changing reference for varIdent parameter estimates in summary.gls

[Kornak, John]

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]]