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
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 observationi
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