Allow for fitting of only fixed effects terms in model formula

[eweine]

Hello,

Thanks for the great software!

I have an application where I would like to fit the following two-group model:

$y_{1}, \dots, y_{n} \sim N(x_{i}^{T} \beta + \alpha T_{i}, \sigma_{I(T_{i} = 1)}^{2} + \sigma_{I(T_{i} = 0)}^{2} + \tau_{i}),$

where $x_{i}$ is a set of covariates, and $T_{i}$ is a binary treatment indicator. The model above has a group specific variance structure, in addition to a known precision component $\tau_{i}$ for each observation.

It seems like your package will allow me to fit the model above when I have a random effect component in the mean formula. However, for my application, I do not.

For example, consider the following simulated data:

library(galamm)

n <- 250
t <- sample(c(0, 1), n, replace = TRUE)
t_var <- ifelse(t == 1, 2, 1)
tau <- rexp(n, rate = 5)
y <- rnorm(n, mean = -0.5 + 0.75 * t, sd = sqrt(t_var + tau))

df <- data.frame(
  y = y,
  t = t,
  tau = tau
)

mod <- galamm(
  y ~ t,
  data = df,
  weights = ~(1|t) + offset(tau)
)

When I attempt to fit this model, I receive the error:

Error: No random effects terms specified in formula

So, would it be possible to allow a model like this to be estimated? If not, are you aware of another package that could fit this?

Thanks!

[osorensen]

Dear @eweine, thanks a lot for reaching out. And sorry, I did not have e-mail notifications on, so did not see your issue until now. I'll take a look tomorrow and get back to you as soon as possible.

[osorensen]

@eweine, I've now taken a look at this, and unfortunately galamm is not set up to estimate this type of model completely without random effects. The nlme package will do the job, however. Below is an example based on your script:

library(nlme)

n <- 250
t <- sample(c(0, 1), n, replace = TRUE)
t_var <- ifelse(t == 1, 2, 1)
tau <- rexp(n, rate = 5)
y <- rnorm(n, mean = -0.5 + 0.75 * t, sd = sqrt(t_var + tau))

df <- data.frame(
  y = y,
  t = t,
  tau = tau
)

mod <- lme(y ~ t, 
           random = ~ 1 | t,
           weights = varIdent(form = ~ 1 | t),
           data = df)
summary(mod)

When running this, the p-value I get for the fixed effect is NA, so there might be an identifiability problem here. However, it seems to give reasonable estimates otherwise.