writing a wrapper for a linear modeling function [MASS::lm.gls()]

368 views Asked by At

The function MASS::lm.gls fits a linear model using generalized least squares, and returns an object of class "lm.gls", but is has no print, summary or other methods.

I could define these simply by hijacking the methods for "lm" objects

print.lm.gls <- function(object, ...) {
     class(object) <- "lm"
     print(object, ...)
}

summary.lm.gls <- function(object, ...) {
     class(object) <- "lm"
     summary(object, ...)
}

Instead, I tried to write a wrapper for lm.gls to add "lm" as another class. (I realize that this might be dangerous, because not all "lm" methods might be valid for GLS.)

Here's what I tried. It doesn't work, as shown in the example below, but I don't understand why not, or how to do this sort of thing more generally.

lm_gls <- function(formula, data, W, subset, na.action, inverse = FALSE, 
    method = "qr", model = FALSE, x = FALSE, y = FALSE, contrasts = NULL, 
    ...) 
{
    result <- MASS::lm.gls(formula, data, W, subset, na.action, 
        inverse = inverse,  method = method, model = model, x = x, y = y, contrasts = contrasts, 
    ...) 
  class(result) <- c(class(result), "lm")
  result
}

Test case:

library(vcd) # needs vcd_1.3-3+ 
data(Punishment, package="vcd")

pun.lor <- loddsratio(Freq ~ memory + attitude + age + education, data = Punishment)
pun.lor.df <- as.data.frame(pun.lor)

library(MASS)

pun.gls <- lm_gls(LOR ~ as.numeric(age) * as.numeric(education), data=pun.lor.df, 
                  W=vcov(pun.lor), inverse=TRUE, x=TRUE, y=TRUE)

This gives the error:

> pun.gls <- lm_gls(LOR ~ as.numeric(age) * as.numeric(education), data=pun.lor.df, W=vcov(pun.lor), inverse=TRUE, x=TRUE, y=TRUE)
Error in xj[i] : invalid subscript type 'closure'
1

There are 1 answers

0
tmpname12345 On BEST ANSWER

The error you are getting is from either the subset or na.action arguments, which are optional but don't have defaults. Hence, when you call the function without specifying them, they are passed as they are from the wrapper, which is type closure. The simplest solution is to pass everything to lm.gls in the dots like so:

lm_gls <- function(...) 
{
  result <- MASS::lm.gls(...) 
  class(result) <- c(class(result), "lm")
  result
}