psychometriko - 10 months ago 41

R Question

**Background:**

I wrote a wrapper function that allows the user to specify the names of variables in a data set to use as predictors in NLME's

`lme()`

`predictors <- str_c(c(w,x), collapse = " + ")`

mod = lme(fixed = reformulate(termlabels = predictors, response = y),

random = reformulate(termlabels = paste0("1|", cluster_label)),

data = dat)

where

`w`

`x`

`cluster_label`

`w = c("w1","w2)`

`x = "x1"`

`school`

`mod = lme(fixed = y ~ w1 + w2 + x1, random = ~1|school, data = tmp)`

The function runs just fine. My problem is that I want to be able to use

`anova()`

`reformulate`

`Error in reformulate(termlabels = predictors, response = y) :`

object 'predictors' not found

and, indeed, when I call

`summary(mod)`

`Fixed: reformulate(termlabels = predictors, response = y)`

I notice that

`anova()`

`...`

These are HSB12 data that reproduce the example just fine. This requires the

`nlme`

`stringr`

`fitVAM <- function(dat,y,w,x,cluster_label) {`

library(nlme)

library(stringr)

predictors <- str_c(c(w,x),collapse=" + ")

mod = lme(fixed = reformulate(termlabels = predictors, response = y),

random = reformulate(termlabels = paste0("1|", cluster_label)),

data = dat)

return(mod)

}

dat <- read.csv("http://www.ats.ucla.edu/stat/paperexamples/singer/hsb12.csv")

mod1 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = "female", cluster_label = "school")

mod2 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = c("female", "ses"), cluster_label = "school")

anova(mod1,mod2)

Answer

There might be a more straight-forward way, but you could use `do.call`

:

```
fitVAM <- function(dat,y,w,x,cluster_label) {
library(nlme)
library(stringr)
predictors <- str_c(c(w,x),collapse=" + ")
params <- list(fixed = reformulate(termlabels = predictors, response = y),
random = reformulate(termlabels = paste0("1|", cluster_label)),
data = dat)
mod = do.call(lme, params)
mod$call[[3]] <- substitute(dat) #otherwise dat is shown expanded in summary and print output
return(mod)
}
#dat <- read.csv("http://www.ats.ucla.edu/stat/paperexamples/singer/hsb12.csv")
mod1 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = "female", cluster_label = "school")
mod2 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = c("female", "ses"), cluster_label = "school")
anova(mod1,mod2)
# Model df AIC BIC logLik Test L.Ratio p-value
#mod1 1 5 46908.59 46942.99 -23449.29
#mod2 2 6 46530.02 46571.29 -23259.01 1 vs 2 380.5734 <.0001
#Warning message:
#In anova.lme(mod1, mod2) :
# fitted objects with different fixed effects. REML comparisons are not meaningful.
```

Source (Stackoverflow)