How to extract weights argument from "gls" object?

Question

I'm using the nlme package to create a generalized least squares model. Here's a reproducible example with a generated dataset:

# Load necessary library
library(nlme)

# Generate some data
set.seed(123)
my_data <- data.frame(
  y = rnorm(100),
  x = runif(100),
  z = factor(rep(1:2, each = 50)),
  g = factor(rep(1:2, each = 50)),
  h = factor(rep(1:2, 50))
)

# Create the model
mdl <- gls(y ~ x * z,
           weights = varIdent(form = ~ 1 | g * h),
           data = my_data)

# I can retrieve the model formula like this
model_formula <- formula(mdl)

In this code, I can use formula(mdl) to get the model formula. However, I can't find a way to retrieve the weights argument from the mdl object. How can I do this?

Answer 1

To retrieve the call:

as.list(mdl$call)$weights
# varIdent(form = ~1 | g * h)

To retrieve the weights you can use the nlme::varWeights() function:

The inverse of the standard deviations corresponding to the variance function structure represented by object are returned.

varWeights(mdl$modelStruct)
#       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2 
# 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 
#       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2 
# 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 
# etc

Answer 2

They're stored as an attribute.

attr(mdl$modelStruct$varStruct, 'weights')
#       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2 
# 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 
#       1*1       1*2       1*1       1*2       1*1       1*2       1*1       1*2 
# 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094 
# ...

str(mdl) is of big help in such cases.

Actually, we could give "gls" a new method.

weights.gls <- function (object, ...) {
  wts <- attr(object$modelStruct$varStruct, 'weights')
  if (is.null(wts)) 
    wts
  else napredict(object$na.action, wts)
}

weights(mdl)
#       1*1       1*2       1*1       1*2       1*1       1*2  ...
# 1.0000000 0.8064094 1.0000000 0.8064094 1.0000000 0.8064094  ...
# ...