Partial R2 in mixed models

Good afternoon;

First, I am not a statistician although I am in the way (I am a medical
doctor studying the grade in statistic, still in the first course). I would
like to compute de partial R-squared of the fixed effects of a model. I
have found a function from the LMERConvenienceFunctions package, but it
computes these from the lme4 anova extraction function, which gives a
sequential anova.

I have created an ad hoc function to compute the R-squared for each term
conditionally to the other terms in the model (based on the pamer.fnc). For
other hand, I have done something similar based with the recommedations
given by Snijders in his book (2nd edition, pages 111-113) to compute de R2
in two levels models.

I am not very sure of what I have done, but I think that the function
works so I would appreciate some light. For other hand, could I call the
calculated value as partial R-squared value of the fixed effects of a mixed
model?

Thank you very much.

As example:

partialR2 <- function(model){
  # Based on SS

  term <- attr(terms(model), "term.labels")
  dv <- gsub(" ", "", gsub("(.*)~.*", "\\1", as.character(model at call)[2]))
  ss.tot <- sum((model at frame[, dv] - mean(model at frame[, dv]))^2)
  n <- length(term)
  ss.var <- numeric(n)
  form <- unlist(lapply(term, function(x) paste(paste(".~. -", x, sep =
""), x, sep = "+")))
  for (i in 1:(n - 1)){
    ss.var[i] <- as.data.frame(anova(update(model,
as.formula(form[i]))))[n, 2]
  }
  ss.var[n] <- as.data.frame(anova(model))[n, 2]
  names(ss.var[n]) <- term[n]
  out <- cbind(round(100 * ss.var/ss.tot, 5))
  rownames(out) <- term
  colnames(out) <- "Partial R2"

  #Snijder
  form <- paste(paste(".~ 1 + (1 |", names(ngrps(model)), sep = ""), ")",
sep = "")
  m.null <- update(model, as.formula(form))
  var.g.null <- VarCorr(m.null)[[1]][1]
  var.r.null <- sigma(m.null)^2
  var.null <- var.g.null + var.r.null

  var.g.full <- VarCorr(model)[[1]][1]
  var.r.full <- sigma(model)^2
  var.full <- var.g.full + var.r.full

  form <- unlist(lapply(term, function(x) paste(".~. -", x, sep = " ")))
  var.red <- numeric(n)
  for (i in n:1){
    var.g.red  <- VarCorr(update(model, as.formula(form[i])))[[1]][1]
    var.r.red <- sigma(update(model, as.formula(form[i])))^2
    var.red[i] <- var.g.red + var.r.red
  }

  out2 <- round(100 * (var.red - var.full)/var.null, 5)
  return(cbind(out, Partial_R2_Snijders = out2))
}


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


library(LMERConvenienceFunctions)
library(foreign)
library(lme4)
dd <- read.dta("
http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_snijders/mlbook1.dta")
str(dd)
m1 <- lmer(langpost ~ sex + ses + iq_perf + langpret + (1|schoolnr), data
= dd)
summary(m1)
anova(m1)
pamer.fnc(m1)
partialR2(m1)

Sorry, yesterday I tried to generalize the function without get it
completely, so I have corrected some mistakes.

Below there is a corrected version.

And, please, apologize my audacity.

=======================================

partialR2 <- function(model){
  # Based on SS
  term <- attr(terms(model), "term.labels")
  dv <- gsub(" ", "", gsub("(.*)~.*", "\\1", as.character(model at call)[2]))
  ss.tot <- sum((model at frame[, dv] - mean(model at frame[, dv]))^2)
  n <- length(term)
  ss.var <- numeric(n)
  form <- unlist(lapply(term, function(x) paste(paste(".~. -", x, sep =
""), x, sep = "+")))

  inter.term <- FALSE

  if (sum(unlist(sapply(term, function(x) grep(":", x)))) >= 1 |
        sum(unlist(sapply(term, function(x) grep("*", x)))) >= 1)
inter.term = TRUE

  if (inter.term == TRUE){
    for (i in 1:(n-1)){
      sub.model <- update(model, as.formula(form[i]))
      ss.var[i] <- as.data.frame(anova(sub.model))[n-1, 2]
    }
    ss.var[n] <- as.data.frame(anova(sub.model))[n, 2]
  }else{
    for (i in 1:(n)){
      sub.model <- update(model, as.formula(form[i]))
      ss.var[i] <- as.data.frame(anova(sub.model))[n, 2]
    }
  }

  names(ss.var[n]) <- term[n]
  out <- cbind(round(100 * ss.var/ss.tot, 5))
  rownames(out) <- term
  colnames(out) <- "Partial R2"

  #Snijder (it gives zero for the interaction components, to find the
values
  #update the models without them.

  form <- paste(paste(".~ 1 + (1 |", names(ngrps(model)), sep = ""), ")",
sep = "")
  m.null <- update(model, as.formula(form))
  var.g.null <- VarCorr(m.null)[[1]][1]
  var.r.null <- sigma(m.null)^2
  var.null <- var.g.null + var.r.null

  var.g.full <- VarCorr(model)[[1]][1]
  var.r.full <- sigma(model)^2
  var.full <- var.g.full + var.r.full

  form <- unlist(lapply(term, function(x) paste(".~. -", x, sep = " ")))
  var.red <- numeric(n)
  for (i in n:1){
    var.g.red  <- VarCorr(update(model, as.formula(form[i])))[[1]][1]
    var.r.red <- sigma(update(model, as.formula(form[i])))^2
    var.red[i] <- var.g.red + var.r.red
  }

  out2 <- round(100 * (var.red - var.full)/var.null, 5)
  return(cbind(out, Partial_R2_Snijders = out2))
}



2014-11-20 17:00 GMT+01:00 Joan Molibo <joanmolibo at gmail.com>:

Good afternoon;

First, I am not a statistician although I am in the way (I am a medical
doctor studying the grade in statistic, still in the first course). I would
like to compute de partial R-squared of the fixed effects of a model. I
have found a function from the LMERConvenienceFunctions package, but it
computes these from the lme4 anova extraction function, which gives a
sequential anova.

I have created an ad hoc function to compute the R-squared for each term
conditionally to the other terms in the model (based on the pamer.fnc). For
other hand, I have done something similar based with the recommedations
given by Snijders in his book (2nd edition, pages 111-113) to compute de R2
in two levels models.

I am not very sure of what I have done, but I think that the function
works so I would appreciate some light. For other hand, could I call the
calculated value as partial R-squared value of the fixed effects of a mixed
model?

Thank you very much.

As example:

partialR2 <- function(model){
  # Based on SS

  term <- attr(terms(model), "term.labels")
  dv <- gsub(" ", "", gsub("(.*)~.*", "\\1", as.character(model at call
)[2]))
  ss.tot <- sum((model at frame[, dv] - mean(model at frame[, dv]))^2)
  n <- length(term)
  ss.var <- numeric(n)
  form <- unlist(lapply(term, function(x) paste(paste(".~. -", x, sep =
""), x, sep = "+")))
  for (i in 1:(n - 1)){
    ss.var[i] <- as.data.frame(anova(update(model,
as.formula(form[i]))))[n, 2]
  }
  ss.var[n] <- as.data.frame(anova(model))[n, 2]
  names(ss.var[n]) <- term[n]
  out <- cbind(round(100 * ss.var/ss.tot, 5))
  rownames(out) <- term
  colnames(out) <- "Partial R2"

  #Snijder
  form <- paste(paste(".~ 1 + (1 |", names(ngrps(model)), sep = ""), ")",
sep = "")
  m.null <- update(model, as.formula(form))
  var.g.null <- VarCorr(m.null)[[1]][1]
  var.r.null <- sigma(m.null)^2
  var.null <- var.g.null + var.r.null

  var.g.full <- VarCorr(model)[[1]][1]
  var.r.full <- sigma(model)^2
  var.full <- var.g.full + var.r.full

  form <- unlist(lapply(term, function(x) paste(".~. -", x, sep = " ")))
  var.red <- numeric(n)
  for (i in n:1){
    var.g.red  <- VarCorr(update(model, as.formula(form[i])))[[1]][1]
    var.r.red <- sigma(update(model, as.formula(form[i])))^2
    var.red[i] <- var.g.red + var.r.red
  }

  out2 <- round(100 * (var.red - var.full)/var.null, 5)
  return(cbind(out, Partial_R2_Snijders = out2))
}


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


library(LMERConvenienceFunctions)
library(foreign)
library(lme4)
dd <- read.dta("
http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_snijders/mlbook1.dta")
str(dd)
m1 <- lmer(langpost ~ sex + ses + iq_perf + langpret + (1|schoolnr), data
= dd)
summary(m1)
anova(m1)
pamer.fnc(m1)
partialR2(m1)