[R-sig-ME] Partial R2 in mixed models

[Joan Molibo]

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

	[[alternative HTML version deleted]]