[R] lme: null deviance, deviance due to the random effects, residual deviance

[Andrew Robinson]

I often want to know what proportion of the variation is being
contributed by different levels of random effects. 

When the random effects are only intercepts, with no slopes, then you
can compute the intra-class correlation, which is the proportion of
the variation explained by the different levels of random effects,
using components from the output of VarCorr.  I wrote a function that
extracts the components and computes them for such models:

varianceDetails <- function(model) {
      if (class(model) != "lme") 
          stop("Input must be an lme!")
      dimension <- floor(dim(VarCorr(model))[1])/2
      j <- c(1:dimension) * 2
      if (length(j) == 1) 
          j <- 1
      i <- c(j, dim(VarCorr(model))[1])
      variances <- as.numeric(VarCorr(model)[i, 1])
      mse <- sum(variances)
      icc <- rep(0, length(j))
      for (k in 1:length(j)) icc[k] <- as.numeric(VarCorr(model)[j[k], 
          1])/mse
      rmse <- sqrt(mse)
      list(icc = icc, rmse = rmse)
}

e.g.

> require(nlme)
> fm1 <- lme(distance ~ age, data = Orthodont)
> varianceDetails(fm1)

(Any feedback on that function is welcomed!)

This function uses the ReML estimates of the variance components (if
the fiting method is ReML!).  When the random effects have slopes then
the problem is more tricky, as the proportion of variance explained
will depend on the predictor variables (see Goldstein et al. at

http://www.mlwin.com/team/materials/pvmm.pdf

).  As a loose rule of thumb, I tend to take a method-of-moments type
of estimate of the standard deviation of the residuals at each level,
and compare those.

Thus,

> sd(Orthodont$distance)
> sd(residuals(fm1, level=0, type="response"))
> sd(residuals(fm1, level=1, type="response"))

You can compute an ad-hoc icc from those components, if you are
comfortable averaging across the random effects.  They will not be the
same as those from varianceDetails, above, because they're not based
on the ReML estimates of the variance components.

I suspect that fitting more complex random effects structures will
render these estimates invalid.  I welcome any further comments on
these ideas.

Cheers

Andrew

> Patrick Giraudoux wrote:
> 
> A maybe trivial and stupid question:
> 
> In the case of a lm or glm fit, it is quite informative (to me) to have 
> a look to the null deviance and the residual deviance of a model. This 
> is generally provided in the print method or the summary, eg:
> 
> Null Deviance:      658.8
> Residual Deviance: 507.3
> 
> and (a bit simpled minded) I like to think that the proportion of 
> deviance 'explained' by the model is (658.8-507.3)/658.8 = 23%
> 
> In the case of lme models, is it possible and reasonable to try and get the:
> - null deviance
> - the total deviance due to the the random effect(s)
> - the residual deviance?
> 
> With the idea that Null deviance = Fixed effects + Random Effects + 
> Residuals
> 
> If yes how to do it ? A lme object provides the following:
> 
>  > names(glm6)
>  [1] "modelStruct"  "dims"         "contrasts"    "coefficients"
>  [5] "varFix"       "sigma"        "apVar"        "logLik"     
>  [9] "numIter"      "groups"       "call"         "method"     
> [13] "fitted"       "residuals"    "fixDF"        "family"
> 
> so no $null.deviance and $deviance elements as in glm objects...
> 
> I tried to find out an answer on R-help & Pineihro & Bates (2000). 
> Partial success only:
> 
> - null deviance: Response: possibly yes: see 
> http://tolstoy.newcastle.edu.au/R/help/05/12/17796.html (Spencer 
> Graves). The (null?) deviance is -2*logLik(mylme), but a personnal trial 
> with some glm objects did not led to the same numbers that the one given 
> by the print.glm method...
> 
> - the deviance due to the the random effect(s). I was supposing that the 
> coefficients given by ranef(mylme) may be an entry... but beyond this, I 
> guess those coefficients must be weighed in some way... which is a far 
> beyond my capacities in this matter...
> 
> - residual deviance. I was supposing that it may be 
> sum(residuals(mylme)^2). With some doubts as far as I feel that I am 
> thinking sum of squares estimation in the context of likelihood and 
> deviance estimations... So  most likely irrelevant. Moreover, in the 
> case I was exploring, this quantity is much larger than the null 
> deviance computed as above...
> 
> Any hint appreciated,
> 
> Patrick Giraudoux
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
Email: a.robinson at ms.unimelb.edu.au         http://www.ms.unimelb.edu.au