Bug 18472 - gls() with singular.ok=TRUE

[hturner]

In the report for Bug 18472 the reporter describes a scenario where the model results in a singular matrix.

Can you create a "reprex", i.e. a reproducible example that has the described behaviour?

In Bug 18285 @bbolker proposed a way to handle rank-deficient lme models and asks if the approach should be applied to other models in the nlme package. If time, you could also look at adapting the approach for gls models - would it handle your reprex more gracefully?

Note that nlme is maintained on https://svn.r-project.org/R-packages/trunk/nlme/ so you would need Subversion to check out the source code for local modification.

[bbolker]

Here's a reprex (I used glmmTMB::simulate_new() for convenience, could re-do without it if including it is a hassle):

library(glmmTMB)
dd <- expand.grid(period = c("before", "after"), ttt = c("control", "treatment"), obs = 1:10)
dd$y <- simulate_new( ~ period*ttt,
                                         dispformula = ~ttt,
                                         family = gaussian,
                                         newdata = dd,
                                         seed = 101,
                                         newparams = list(beta = c(1,2,1,2), betad = c(0,1)))[[1]]
dd_lim <- subset(dd, !(period == "before" & ttt == "treatment"))
library(nlme)
fit <- gls(y ~ period*ttt, weights = varIdent(form = ~1|ttt), data  = dd)
fit2 <- update(fit, data = dd_lim,  control = glsControl(singular.ok = TRUE))

summary() applied to both shows that the weights are not estimated for fit2.

[hturner]

Thanks Ben. This serves the purpose for initial testing but it would be better to have an example that doesn't rely on additional packages to share on Bugzilla.

[bbolker]

OK. I think we could replace simulate_new with:

set.seed(101)
X <- model.matrix(~period*ttt, data = dd)
mu <- X %*% c(1,2,1,2)
sd <- ifelse(dd$ttt == "control", 1, exp(1))
dd$y <- rnorm(nrow(dd), mean = mu, sd = sd)

(not tested).

[hturner]

Bug 18472 was updated with the reprex that could be fixed using the approach Ben proposed in Bug 18285. However it raised issues with the proposed interface, so a revised proposal was made on Bug 18285.