rdiagnosticsnlme

Calculqting leverage and Cook's distance in non-linear mixed effects model


I am trying to calculate leverage and cooks distance for the model that I developed using nlme, but I'm getting the error saying

Error in cooks.distance.lme(model) : not implemented for "nlme" objects"

What are the ways to calculate Cook's distance in nlme?

formula = log_Lobs ~ log(150*((1 + ((150/Lt_1)^(1/exp(p))-1)*exp(-exp(k)*td/365))^(-exp(p))))
model <- do.call(nlme,
                    list(formula,
                         fixed = c(p ~ 1, k ~ 1 + season2),
                         random = k ~ 1 | id,
                         data = data_select,
                         start = list(fixed = c(p, k)),
                         na.action = na.exclude,
                         control=list(maxIter=1e6, msMaxIter = 1e6, msVerbose = TRUE)
                    ))

#Calculate leverage
lev<-hat(model.matrix(model))

Error in x$terms %||% attr(x, "terms") %||% stop("no terms component nor attribute") : no terms component nor attribute

#Calculate Cook's Distance
cd<-cooks.distance(model)

Error in cooks.distance.lme(model) : not implemented for "nlme" objects


Solution

  • You could do it by brute force: from Wikipedia,

    Cook's distance of observation i is defined as the sum of all the changes in the regression model when observation i is removed from it

    Using the formula from Nieuwenhuis 2012

    library(nlme)
    mod1 <- lme(distance ~ age*Sex, random = ~age|Subject, data = Orthodont)
    coef1 <- fixef(mod1)
    cooksd <- rep(NA, nrow(Orthodont))
    Vinv <- solve(vcov(mod1))
    np <- length(coef1)
    for (i in seq(nrow(Orthodont))) {
       newmod <- try(update(mod1, data = Orthodont[-i, ]), silent = TRUE)
       if (inherits(newmod, "try-error")) next
       newcoef <- fixef(newmod)
       cdiff <- coef1-newcoef
       cooksd[i] <- (cdiff %*% Vinv %*% cdiff)/np
    }
    

    Nieuwenhuis, Rense. 2012. “Influence.ME: Tools for Detecting Influential Data in Mixed Effects Models” R Journal 4/2 https://journal.r-project.org/archive/2012/RJ-2012-011/RJ-2012-011.pdf