Search code examples
rregressionlmmixed-modelsnlme

Shortening the formula syntax of a regression model


I was wondering if the syntax of the regression model below could be made more concise (shorter) than it currently is?

dat <- read.csv('https://raw.githubusercontent.com/rnorouzian/v/main/bv1.csv')

library(nlme)

model <- lme(achieve ~ 0 + D1 + D2+ 
       D1:time   + D2:time+
       D1:schcontext + D2:schcontext + 
       D1:female + D2:female+
       D1:I(female*time) + D2:I(female*time)+
       D1:I(schcontext*time) + D2:I(schcontext*time), correlation = corSymm(),
       random = ~0 + D1:time | schcode/id, data = dat, weights = varIdent(form = ~1|factor(math)),
     na.action = na.omit, control = lmeControl(maxIter = 200, msMaxIter = 200, niterEM = 50,
                                               msMaxEval = 400))
coef(summary(model))

Solution

  • Focusing on the fixed-effect component only.

    Original formula:

    form1 <-  ~ 0 + D1 + D2+ 
           D1:time   + D2:time+
           D1:schcontext + D2:schcontext + 
           D1:female + D2:female+
           D1:I(female*time) + D2:I(female*time)+
           D1:I(schcontext*time) + D2:I(schcontext*time)
    X1 <- model.matrix(form1, data=dat) 
    

    I think this is equivalent

    form2 <- ~0 +
        D1 + D2 + 
        (D1+D2):(time + schcontext + female + female:time+schcontext:time)
    X2 <- model.matrix(form2, data=dat)
    

    (Unfortunately ~ 0 + (D1 + D2):(1 + time + ...) doesn't work as I would have liked/expected.) For a start, the model matrix has the right dimensions. Staring at the column names of the model matrices and reordering the columns manually:

    X2o <- X2[,c(1:3,6,4,7,5,8,9,11,10,12)]
    all.equal(c(X1),c(X2o)) ##TRUE
    

    (For numerical predictors, you don't need I(A*B): A:B is equivalent.)

    Actually you can do a little better using the * operator

    form3 <- ~0 +
        D1 + D2 + 
        (D1+D2):(time*(schcontext+female))
    X3 <- model.matrix(form3, data=dat)
    X3o <- X3[,c(1:3,6,4,7,5,8,10,12,9,11)]
    all.equal(c(X1),c(X3o))  ## TRUE
    

    Compare formula length:

    sapply(list(form1,form2,form3),
           function(x) nchar(as.character(x)[[2]]))
    ## [1] 183 84 54