Implement a refit function

[florianhartig]

Hi Dimitris,

I am currently making another attempt to integrate GLMMadaptive to DHARMa, see https://github.com/florianhartig/DHARMa/issues/96

One point that I noted is that GLMMadaptive does not have a refit function, as in https://www.rdocumentation.org/packages/lme4/versions/1.1-26/topics/refit

I tried to circumvent this via the following code

getRefit.MixMod <- function(object, newresp, ...) { responsename = colnames(model.frame(object))[1] newDat = object$data newDat[, match(responsename,names(newDat))] = newresp update(object, data = newDat) }

The issue with that is that this fails for k/n binomial with c(s,f) ~ pred syntax

Any great idea how to solve this with a quick hack that is not super hacky, or do you see an option to supply a refit function in GLMMadaptive?

[drizopoulos]

Hi Florian,

Thanks for taking this up. I will have a look and let you know how this could be done.

However, it is not clear to me why you need to refit the model. This will be computationally expensive. My understanding was that you only need to simulate from the model. At least, this is how I use your package in the custom made function resids_plot() I have here: https://drizopoulos.github.io/GLMMadaptive/articles/Goodness_of_Fit.html.

Best, Dimitris

From: Florian Hartig notifications@github.com Sent: Monday, January 25, 2021 10:54 PM To: drizopoulos/GLMMadaptive GLMMadaptive@noreply.github.com Cc: Subscribed subscribed@noreply.github.com Subject: [drizopoulos/GLMMadaptive] Implement a refit function (#32)

Hi Dimitris,

I am currently making another attempt to integrate GLMMadaptive to DHARMa, see florianhartig/DHARMa#96https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fflorianhartig%2FDHARMa%2Fissues%2F96&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cfabf1c34d9e54aeb9f7408d8c17bc033%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637472084536138576%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=TrVLcX3eWCmgzd3KCkNwb9F3HjFWggDNKMseW9aRA5A%3D&reserved=0

One point that I noted is that GLMMadaptive does not have a refit function, as in https://www.rdocumentation.org/packages/lme4/versions/1.1-26/topics/refithttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.rdocumentation.org%2Fpackages%2Flme4%2Fversions%2F1.1-26%2Ftopics%2Frefit&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cfabf1c34d9e54aeb9f7408d8c17bc033%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637472084536138576%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=F7VlZk86BPK5wwofbIpLAjyfS50IUlPEMpVDUoVdlaE%3D&reserved=0

I tried to circumvent this via the following code

getRefit.MixMod <- function(object, newresp, ...) { responsename = colnames(model.frame(object))[1] newDat = object$data newDat[, match(responsename,names(newDat))] = newresp update(object, data = newDat) }

The issue with that is that this fails for k/n binomial with c(s,f) ~ pred syntax

Any great idea how to solve this with a quick hack that is not super hacky, or do you see an option to supply a refit function in GLMMadaptive?

- You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrizopoulos%2FGLMMadaptive%2Fissues%2F32&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cfabf1c34d9e54aeb9f7408d8c17bc033%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637472084536138576%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=z%2BmWNgHYHB9OM391Hho5wqozyxMi%2BAJsBYfwNT%2B1Ta8%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADE7TTYHJKS5PPIY2TXMSF3S3XSAFANCNFSM4WSNKVFQ&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cfabf1c34d9e54aeb9f7408d8c17bc033%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637472084536148534%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=0BOINlFKbM4nGZJVvwR2EvlJI7tabajG4qmW%2F6n5QxM%3D&reserved=0.

[florianhartig]

Hi Dimitris,

for historical reasons, I'm still carrying around 2 versions of calculating everything in DHARMa,

1) where quantile residuals are created by comparing observed vs. simulated data (default) 2) where the models are re-fit to the simulated data, and residuals or other properties of the model are compared

The original motivation for this was the idea that option 2 would be more reliable for biased estimators. In practice, I have found few advantages of the refit, but for the moment I'm keeping the option.

More generally though, the refit is important for parametric bootstrap and all that.

[drizopoulos]

Hi Florian,

I have written a short refit() function that seems to do the trick - see below:

# simulate some data
set.seed(123L)
n <- 500
K <- 15
t.max <- 25

betas <- c(-2.13, -0.25, 0.24, -0.05)
D <- matrix(0, 2, 2)
D[1:2, 1:2] <- c(0.48, -0.08, -0.08, 0.18)

times <- c(replicate(n, c(0, sort(runif(K-1, 0, t.max)))))
group <- sample(rep(0:1, each = n/2))
DF <- data.frame(year = times, group = factor(rep(group, each = K)))
X <- model.matrix(~ group * year, data = DF)
Z <- model.matrix(~ year, data = DF)

b <- cbind(rnorm(n, sd = sqrt(D[1, 1])), rnorm(n, sd = sqrt(D[2, 2])))
id <- rep(1:n, each = K)
eta.y <- as.vector(X %*% betas + rowSums(Z * b[id, ]))
DF$y <- rbinom(n * K, 10, plogis(eta.y))
DF$fails <- 10 - DF$y
DF$id <- factor(id)

################################################################
################################################################

fm <- mixed_model(fixed = cbind(y, fails) ~ year * group, 
                  random = ~ 1 | id, data = DF,
                  family = binomial())

refit_MixMod <- function (object, newresp = NULL) {
    if (is.null(newresp)) return(object)
    if (!is.data.frame(newresp)) stop("'newresp' must be a data.frame.")
    newdata <- object$data
    respVars <- all.vars(update.formula(formula(object), . ~ 0))
    newdata[respVars] <- newresp[respVars]
    update(object, data = newdata)
}

fm
refit_MixMod(fm, DF)