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

[Joan Molibo]

Finally, I have found this work An R^2 statistics for fixed effects in the
linear mixed model. Stat Med 2008 28: 6137, which has been of great help.

R2.mmixed <- function(model, d = 5){
  require(pbkrtest)
  term <- attr(terms(model), "term.labels")
  n.term <- length(term)
  form <- paste(paste(".~ 1 + (1 |", names(ngrps(model)), sep = ""), ")",
sep = "")
  m.null <- update(model, as.formula(form))
  kr <- KRmodcomp(model, m.null)[[1]]
  R2 <- kr[1, 1] * kr[1, 3]^(-1) * length(term)/(1 + kr[1, 1] * kr[1,
3]^(-1) * n.term)
  R2 <- round(R2 * 100, d)

  aov.model <- Anova(model, test.statistic = "F", type = 3)
  p.R2 <- numeric(n.term)
  for (i in 1:n.term){
    p.R2[i] <- aov.model[i+1, 1] * (aov.model[i+1, 3])^(-1) *
aov.model[i+1, 2]/(1 + aov.model[i+1, 1] * (aov.model[i+1, 3])^(-1) *
aov.model[i+1, 2])
  }
  p.R2 <- matrix(round(p.R2*100, d), ncol = 1);
  rownames(p.R2) <- term
  colnames(p.R2) <- "Estimate"
  l <- list("Model R-squared (fix effects)" = R2, "Partial R-squared" =
p.R2)
  return(l)
}

2014-11-21 12:06 GMT+01:00 Joan Molibo <joanmolibo at gmail.com>:

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