Search code examples
rinteractionmixed-models

R: cant get a lme{nlme} to fit when using self-constructed interaction variables

I'm trying to get a lme with self constructed interaction variables to fit. I need those for post-hoc analysis.

library(nlme)

# construct fake dataset
obsr <- 100
dist <- rep(rnorm(36), times=obsr)
meth <- dist+rnorm(length(dist), mean=0, sd=0.5); rm(dist)
meth <- meth/dist(range(meth)); meth <- meth-min(meth)
main <- data.frame(meth = meth,
                   cpgl = as.factor(rep(1:36, times=obsr)),
                   pbid = as.factor(rep(1:obsr, each=36)),
                   agem = rep(rnorm(obsr, mean=30, sd=10), each=36),
                   trma = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)),
                   depr = as.factor(rep(sample(c(TRUE, FALSE), size=obsr, replace=TRUE), each=36)))

# check if all factor combinations are present
# TRUE for my real dataset; Naturally TRUE for the fake dataset
with(main, all(table(depr, trma, cpgl) >= 1))

# construct interaction variables
main$depr_trma <- interaction(main$depr, main$trma, sep=":", drop=TRUE)
main$depr_cpgl <- interaction(main$depr, main$cpgl, sep=":", drop=TRUE)
main$trma_cpgl <- interaction(main$trma, main$cpgl, sep=":", drop=TRUE)
main$depr_trma_cpgl <- interaction(main$depr, main$trma, main$cpgl, sep=":", drop=TRUE)

# model WITHOUT preconstructed interaction variables
form1 <- list(fixd = meth ~ agem + depr + trma + depr*trma + cpgl +
                            depr*cpgl +trma*cpgl + depr*trma*cpgl,
              rndm = ~ 1 | pbid,
              corr = ~ cpgl | pbid)

modl1 <- nlme::lme(fixed=form1[["fixd"]],
                   random=form1[["rndm"]],
                   correlation=corCompSymm(form=form1[["corr"]]),
                   data=main)

# model WITH preconstructed interaction variables
form2 <- list(fixd = meth ~ agem + depr + trma + depr_trma + cpgl +
                            depr_cpgl + trma_cpgl + depr_trma_cpgl,
              rndm = ~ 1 | pbid,
              corr = ~ cpgl | pbid)

modl2 <- nlme::lme(fixed=form2[["fixd"]],
                   random=form2[["rndm"]],
                   correlation=corCompSymm(form=form2[["corr"]]),
                   data=main)

The first model fits without any problems whereas the second model gives me following error:

Error in MEEM(object, conLin, control$niterEM) : 
  Singularity in backsolve at level 0, block 1

Nothing i found out about this error so far helped me to solve the problem. However the solution is probably pretty easy.

Can someone help me? Thanks in advance!

EDIT 1:

When i run:

modl3 <- lm(form1[["fixd"]], data=main)
modl4 <- lm(form2[["fixd"]], data=main)

The summaries reveal that modl4 (with the self constructed interaction variables) in contrast to modl3 shows many more predictors. All those that are in 4 but not in 3 show NA as coefficients. The problem therefore definitely lies within the way i create the interaction variables...

EDIT 2:

In the meantime I created the interaction variables "by hand" (mainly paste() and grepl()) - It seems to work now. However I would still be interested in how i could have realized it by using the interaction() function.

Solution

  • I should have only constructed the largest of the interaction variables (combining all 3 simple variables).

    If i do so the model gets fit. The likelihoods then are very close to each other and the number of coefficients matches exactly.