nlme How to extract and understand the coefficients of the variance function

Question

I am using nlme to fit heteroskedastic residual variance as a function of covariates. I am not sure how to extract and understand the values produced by nlme. Particularly the output of the coef function on the variance function.

Here is some context:

In contexts of experiments where the within residual variance differs in the treatment condition, we can model the heteroskedastic directly as a function of the covariate treatment. Seen generally as a variance function: $$ \text{Var}(e_i|{\bf b}_i) = \sigma^2g^2(\mu_{ij},\boldsymbol{\nu}_{ij}, \boldsymbol{\delta}) $$ Where $\boldsymbol{\nu}_{ij}$ is the variance covariates and $\boldsymbol{\delta}$ are variance parameters.

In my case this simplifies as a different variance in the treatment groups which make up strata $s$

$$ \text{Var}(e_i) = \sigma^2\delta_s^2 $$ In my case treatment is binary so there are 2 strata. Since there are 3 variance parameters with this specification $\sigma$, $\delta_c$ for the control and $\delta_t$ for the treatment. One of the $\delta$ is redundant and set to 1. $\delta$ then represents that ratio of variance between the $\delta_s$ that was set to 1.

In this simulated example where the treatment as double the variance in control group

$$ Y_i = 3 + e^c_i + T_ie^t_i $$ $$ e^c_i,e^t_i \sim N(0,1) $$ $$ T_i \sim Bern(.5) $$

Simulating data. Note the double variance in the treatment group.

library(nlme) #For modeling
library(dplyr) #For data maniplulation
set.seed(123)
n <- 1000
treatment <- rbinom(n = n, size = 1, prob = .5)
y <-  3 + rnorm(n) + rnorm(n)*treatment #treatment will have double variance

data.frame(y, treatment) |> 
  group_by(treatment) |> 
  summarise(var(y))
#> # A tibble: 2 × 2
#>   treatment `var(y)`
#>       <int>    <dbl>
#> 1         0     1.03
#> 2         1     2.13

Fitting heteroskedastic residuals using nlme:

model <- gls(y ~ 1, weights = varIdent(form = ~1|treatment))

model$modelStruct$varStruct #How to extract the values?
#> Variance function structure of class varIdent representing
#>        0        1 
#> 1.000000 1.437717

1.43 represents the value of $\delta_t$ and when squared equals around 2 which is the relationship between the control variance and treatment variance

I am not sure how to extract the values of $\delta$ from the gls object.

Using the coef function on the variance produces a value that I do not know what it is.

coef(model$modelStruct$varStruct) # what does this value mean?
#> [1] 0.3630564

References

Pinheiro, J., & Bates, D. (2006). Mixed-effects models in S and S-PLUS. Springer science & business media. Chapter 5.2

Answer 1

I looked at the print function for objects of variance functions that https://svn.r-project.org/R-packages/trunk/nlme/R/varFunc.R.

Found you need to specify unconstrained = FALSE to get the values that represent $\delta$. When unconstrained = TRUE it returns the values natural log.

model$modelStruct$varStruct
#> Variance function structure of class varIdent representing
#>        0        1 
#> 1.000000 1.437717

coef(model$modelStruct$varStruct, unconstrained = FALSE, allCoef = TRUE)
#>        0        1 
#> 1.000000 1.437717

From the documentation this is the "natural", generally constrained form.

coef(model$modelStruct$varStruct, unconstrained = TRUE, allCoef = TRUE)
#> [1] 0.3630564
exp(coef(model$modelStruct$varStruct, unconstrained = TRUE, allCoef = TRUE))
#> [1] 1.437717

This value is the $\text{ln}(\delta)$