5
$\begingroup$

I am modelling data exhibiting a tweedie distribution in R using glmer (package lme4). To compare the models I would like to use the AIC and R^2. I have a couple of question on this (example code at the end):

1) Should I even be using the tweedie Family with lme4 and glmer?

2) The AICtweedie function (package tweedie) accepts only glm objects, maybe as many don't agree with AIC for glmm (?), but I adapted it to accept a glmer model object - So my question is, assuming AIC is a good idea for glmm would I simply calculate it in the same way as for a glm?

Adapted AICtweedie function:

AICtweedie_glmer <- function(glmer.obj, k = 2){
  wt <- weights(glmer.obj)
  n <- length(residuals(glmer.obj))
  edf <- n - df.residual(glmer.obj)
  dev <- deviance(glmer.obj)
  disp <- dev/n
  mu <- glmer.obj@resp$mu
  y <- glmer.obj@resp$y
  p <- get("p", envir = environment(glmer.obj@resp$family$variance))
  den <- dtweedie(y = y, mu = mu, phi = disp, power = p)
  AIC <- -2 * sum(log(den) * wt)
  return(AIC + k * (edf + 1))
}

3) I wanted to use the r.squaredGLMM function (package piecewiseSEM) to calculate marginal and conditional R^2 but it seems the tweedie Family is not yet implemented. apologies for not inline formulae, I couldn't embed Images.

In the paper of Nakagawa et al. ("The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded") they say the formulations (marginal & conditional) are applicable to tweedie distributions too. The formulae contain the Terms: variance of fixed effects, variance of random effects, and residual variance or observation-level variance .

From the publication I can obtain observation-level variance using this, with μ "the mean on the observed scale" and φ "the dispersion Parameter"

I think observation-level variance can therefore be calculated from a glmer object so:

m <- mean(glmer.obj@resp$y)
p <- get("p", envir = environment(glmer.obj@resp$family$variance))
disp <- dev/n
varObs <- (disp*(m^p)^-2)

Taking the forumla inside the r.squaredGLMM function (marginal for example):

Rm <- varF/(varF+varRand+varDist+varDisp)

So my second question is: Is it correct to simply replace varDist + varDisp with my calculated varObs? Or have I overlooked something?

# EXAMPLE #

# create a glmer object with Response of tweedie dist:
y1 <- data.frame(y= rtweedie(500, p= 1.6, mu= 150, phi= 0.65), x="1")
y2 <- data.frame(y=rtweedie(500, p= 1.6, mu= 200, phi= 0.65), x="2")

glmer.obj <- glmer(y ~ 1 +(1|x),data=rbind(y1,y2),family=tweedie(var.power = 1.6))

# varF and varRand taken from package piecewiseSEM:
varF <- var(as.vector(lme4::fixef(glmer.obj) %*% t(glmer.obj@pp$X)))

varRand <- sum(sapply(
  VarCorr(glmer.obj)[!sapply(unique(unlist(strsplit(names(ranef    (glmer.obj)),":|/"))), 
                  function(l) length(unique(glmer.obj@frame[,l])) == nrow(glmer.obj@frame))],
function(Sigma){X <- model.matrix(glmer.obj)
              Z <- X[,rownames(Sigma)]
              sum(diag(Z %*% Sigma %*% t(Z)))/nrow(X)})) 

n <- length(residuals(glmer.obj))
dev <- deviance(glmer.obj)
m <- mean(glmer.obj@resp$y)
p <- get("p", envir = environment(glmer.obj@resp$family$variance))
disp <- dev/n
varObs <- (disp*(m^p)^-2)


# Calculate marginal R-squared 
Rm <- varF/(varF+varRand+varObs) 

# Calculate conditional R-squared (fixed effects+random effects/total variance) 
Rc <- (varF+varRand)/(varF+varRand+varObs)

Thanks stackoverflow community! (from a long time reader, first time poster).

$\endgroup$

closed as too broad by mdewey, kjetil b halvorsen, Robert Long, gung Aug 19 at 15:14

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

migrated from stackoverflow.com Jan 29 '18 at 16:11

This question came from our site for professional and enthusiast programmers.

  • 1
    $\begingroup$ "I am modelling data exhibiting a tweedie distribution" The distribution of data is not very relevant for regression analysis. What is important is the distribution of residuals. There is currently no way to assess from your question if use of tweedie distribution is appropriate/necessary at all, apart from the implementation details. $\endgroup$ – Roland Jan 29 '18 at 14:41
  • 1
    $\begingroup$ @Roland Thanks for your Response! I thought I should employ tweedie family because my response is Zero inflated - Looks somewhat like the top-right graph here. The residuals after fitting my base model are quite even (I can't post Images as I'm too new), the variance can be assumed to be constant. $\endgroup$ – Joe Pizzle Jan 29 '18 at 15:03