Re: [R] Effect size in mixed models

[Andrew Robinson]

Hi Spencer,

what you're describing now seems to me to be similar to the idea of
intra-class correlation.  I have previously used the relative and
absolute amounts of variance represented at different hierarchical
levels of random effects, and they seemed to me to be easy to
interpret and useful, although I used method of moments estimators
instead of ML.  This ease of interpretation and utility is especially
true when the random effects structure matches the physical layout of
a sample, e.g., multiple measurements within trees within plots within
stands within forests etc.
 
Sadly I cannot see how to cleanly generalize the principle to
situations of crossed random effects.  Also, as Goldstein et al point
out in a nice paper in Understanding Statistics (2002): "Partitioning
variation in multilevel models", random covariates also complicate the
picture substantially.

Cheers

Andrew

On Thu, Jun 22, 2006 at 05:23:52PM -0700, Spencer Graves wrote:
> I just learned that my earlier suggestion was wrong.  It's better to 
> compute the variance of the predicted or fitted values and compare those 
> with the estimated variance components.
> 
> To see how to do this, consider the following minor modification of 
> an example in the "lme" documentation:
> 
> fm1. <- lme(distance ~ age, data = Orthodont, random=~1)
> fm2. <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
> 
> # str(fm1.) suggested the following:
>  > var(fm2.$fitted[, "fixed"]-fm1.$fitted[, "fixed"])
> [1] 1.312756
>  > VarCorr(fm1.)[, 1]
> (Intercept)Residual
>   "4.472056"  "2.049456"
>  > VarCorr(fm2.)[, 1]
> (Intercept)Residual
>   "3.266784"  "2.049456"
> 
> In this example, the subject variance without considering "Sex" was 
> 4.47 but with "Sex" in the model, it dropped to 3.27, while the Residual 
> variance remained unchanged at 2.05.  The difference between 
> fm2.$fitted[, "fixed"] and fm1.$fitted[, "fixed"] is the change in the 
> predictions generated by the addition of "Sex" to the model.  The 
> variance of that difference was 1.31.  Note that 3.27 + 1.31 = 4.58, 
> which is moderately close to 4.47.
> 
> In sum, I think we can get a reasonable estimate of the size of an 
> effect from the variance of the differences in the "fixed" portion of 
> the fitted model.
> 
> Comments?
> Hope this helps.
> Spencer Graves
> 
> Spencer Graves wrote:
> >   You have asked a great question:  It would indeed be useful to 
> > compare the relative magnitude of fixed and random effects, e.g. to 
> > prioritize efforts to better understand and possibly manage processes 
> > being studied.  I will offer some thoughts on this, and I hope if there 
> > are errors in my logic or if someone else has a better idea, we will 
> > both benefit from their comments.
> > 
> >   The ideal might be an estimate of something like a mean square for 
> > a particular effect to compare with an estimated variance component.
> > Such mean squares were a mandatory component of any analysis of variance 
> > table prior to the (a) popularization of generalized linear models and 
> > (b) availability of software that made it feasible to compute maximum 
> > likelihood estimates routinely for unbalanced, mixed-effects models. 
> > However, anova(lme(...)) such mean squares are for most purposes 
> > unnecessary cluster in a modern anova table.
> > 
> >   To estimate a mean square for a fixed effect, consider the 
> > following log(likelihood) for a mixed-effects model:
> > 
> >   lglk = (-0.5)*(n*log(2*pi*var.e)-log(det(W)) + 
> > t(y-X%*%b)%*%W%*%(y-X%*%b)/var.e),
> > 
> > where n = the number of observations,
> > 
> >   b = the fixed-effect parameter variance,
> > 
> > and the covariance matrix of the residuals, after integrating out the 
> > random effects is var.e*solve(W).  In this formulation, the matrix "W" 
> > is a function of the variance components.  Since they are not needed to 
> > compute the desired mean squares, they are suppressed in the notation here.
> > 
> >   Then, the maximum likelihood estimate of
> > 
> >   var.e = SSR/n,
> > 
> > where SSR = t(y-X%*%b)%*%W%*%(y-X%*%b).
> > 
> >   Then
> > 
> >   mle.lglk = (-0.5)*(n*(log(2*pi*SSR/n)-1)-log(det(W))).
> > 
> >   Now let
> > 
> >   SSR0 = this generalized sum of squares of residuals (SSR) without 
> > effect "1",
> > 
> > and
> > 
> >   SSR1 = this generalized SSR with this effect "1".
> > 
> >   If I've done my math correctly, then
> > 
> >   D = deviance = 2*log(likelihood ratio)
> > = (n*log(SSR0/SSR1)+log(det(W1)/det(W0)))
> > 
> >   For roughly half a century, a major part of "the analysis of 
> > variance" was the Pythagorean idea that the sum of squares under H0 was 
> > the sum of squares under H1 plus the sum of squares for effect "1":
> > 
> >   SSR0 = SS1 + SSR1.
> > 
> >   Whence,
> > 
> >   exp((D/n)-log(det(W1)/det(W0))) = 1+SS1/SSR1.
> > 

Re: [R] Effect size in mixed models

[Dave Atkins]

Spencer & Bruno--

You might to take a look at a paper by Kyle Roberts at:

http://www.hlm-online.com/papers/

Kyle has been working on the issue of effect-sizes in mixed-effects (aka 
multilevel aka HLM) for a couple years.  I haven't had a chance to compare what 
Spencer has suggested with Kyle's approach.  [Though, it would be *lovely* if 
there were an agreed upon method for estimating effect-sizes in mixed-effects 
models; in psychology, reviewers will bludgeon you for an effect-size...]

cheers, Dave
-- 
Dave Atkins, PhD
Assistant Professor in Clinical Psychology
Fuller Graduate School of Psychology
Email: [EMAIL PROTECTED]



Spencer wrote:

   I just learned that my earlier suggestion was wrong.  It's better to
compute the variance of the predicted or fitted values and compare those
with the estimated variance components.

  To see how to do this, consider the following minor modification of
an example in the "lme" documentation:

fm1. <- lme(distance ~ age, data = Orthodont, random=~1)
fm2. <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)

# str(fm1.) suggested the following:
  > var(fm2.$fitted[, "fixed"]-fm1.$fitted[, "fixed"])
[1] 1.312756
  > VarCorr(fm1.)[, 1]
(Intercept)Residual
   "4.472056"  "2.049456"
  > VarCorr(fm2.)[, 1]
(Intercept)Residual
   "3.266784"  "2.049456"

  In this example, the subject variance without considering "Sex" was
4.47 but with "Sex" in the model, it dropped to 3.27, while the Residual
variance remained unchanged at 2.05.  The difference between
fm2.$fitted[, "fixed"] and fm1.$fitted[, "fixed"] is the change in the
predictions generated by the addition of "Sex" to the model.  The
variance of that difference was 1.31.  Note that 3.27 + 1.31 = 4.58,
which is moderately close to 4.47.

  In sum, I think we can get a reasonable estimate of the size of an
effect from the variance of the differences in the "fixed" portion of
the fitted model.

  Comments?
  Hope this helps.
  Spencer Graves

Spencer Graves wrote:
 >   You have asked a great question:  It would indeed be useful to
 > compare the relative magnitude of fixed and random effects, e.g. to
 > prioritize efforts to better understand and possibly manage processes
 > being studied.  I will offer some thoughts on this, and I hope if there
 > are errors in my logic or if someone else has a better idea, we will
 > both benefit from their comments.
 >
 >   The ideal might be an estimate of something like a mean square for
 > a particular effect to compare with an estimated variance component.
 > Such mean squares were a mandatory component of any analysis of variance
 > table prior to the (a) popularization of generalized linear models and
 > (b) availability of software that made it feasible to compute maximum
 > likelihood estimates routinely for unbalanced, mixed-effects models.
 > However, anova(lme(...)) such mean squares are for most purposes
 > unnecessary cluster in a modern anova table.
 >
 >   To estimate a mean square for a fixed effect, consider the
 > following log(likelihood) for a mixed-effects model:
 >
 >   lglk = (-0.5)*(n*log(2*pi*var.e)-log(det(W)) +
 > t(y-X%*%b)%*%W%*%(y-X%*%b)/var.e),
 >
 > where n = the number of observations,
 >
 >   b = the fixed-effect parameter variance,
 >
 > and the covariance matrix of the residuals, after integrating out the
 > random effects is var.e*solve(W).  In this formulation, the matrix "W"
 > is a function of the variance components.  Since they are not needed to
 > compute the desired mean squares, they are suppressed in the notation here.
 >
 >   Then, the maximum likelihood estimate of
 >
 >   var.e = SSR/n,
 >
 > where SSR = t(y-X%*%b)%*%W%*%(y-X%*%b).
 >
 >   Then
 >
 >   mle.lglk = (-0.5)*(n*(log(2*pi*SSR/n)-1)-log(det(W))).
 >
 >   Now let
 >
 >   SSR0 = this generalized sum of squares of residuals (SSR) without
 > effect "1",
 >
 > and
 >
 >   SSR1 = this generalized SSR with this effect "1".
 >
 >   If I've done my math correctly, then
 >
 >   D = deviance = 2*log(likelihood ratio)
 > = (n*log(SSR0/SSR1)+log(det(W1)/det(W0)))
 >
 >   For roughly half a century, a major part of "the analysis of
 > variance" was the Pythagorean idea that the sum of squares under H0 was
 > the sum of squares under H1 plus the sum of squares for effect "1":
 >
 >   SSR0 = SS1 + SSR1.
 >
 >   Whence,
 >
 >   exp((D/n)-log(det(W1)/det(W0))) = 1+SS1/SSR1.
 >
 > Thus,
 >
 >   SS1 = SSR1*(exp((D/n)-log(det(W1)/det(W0)))-1).
 >
 >   If the difference between deg(W1) and det(W0) can be ignored, we get:
 >
 >   SS1 = SSR1*(exp((D/n)-1).
 >
 >   Now compute MS1 = SS1/df1, and compare with the variance components.
 >
 >   If there is a flaw in this logic, I hope someone will disabuse me
 > of it.
 >
 >   If this seems too terse or convoluted to follow, please provide a
 > simple, self-contained example, as sugges

Re: [R] Effect size in mixed models

[Spencer Graves]

  I just learned that my earlier suggestion was wrong.  It's better to 
compute the variance of the predicted or fitted values and compare those 
with the estimated variance components.

  To see how to do this, consider the following minor modification of 
an example in the "lme" documentation:

fm1. <- lme(distance ~ age, data = Orthodont, random=~1)
fm2. <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)

# str(fm1.) suggested the following:
 > var(fm2.$fitted[, "fixed"]-fm1.$fitted[, "fixed"])
[1] 1.312756
 > VarCorr(fm1.)[, 1]
(Intercept)Residual
  "4.472056"  "2.049456"
 > VarCorr(fm2.)[, 1]
(Intercept)Residual
  "3.266784"  "2.049456"

  In this example, the subject variance without considering "Sex" was 
4.47 but with "Sex" in the model, it dropped to 3.27, while the Residual 
variance remained unchanged at 2.05.  The difference between 
fm2.$fitted[, "fixed"] and fm1.$fitted[, "fixed"] is the change in the 
predictions generated by the addition of "Sex" to the model.  The 
variance of that difference was 1.31.  Note that 3.27 + 1.31 = 4.58, 
which is moderately close to 4.47.

  In sum, I think we can get a reasonable estimate of the size of an 
effect from the variance of the differences in the "fixed" portion of 
the fitted model.

  Comments?
  Hope this helps.
  Spencer Graves

Spencer Graves wrote:
>   You have asked a great question:  It would indeed be useful to 
> compare the relative magnitude of fixed and random effects, e.g. to 
> prioritize efforts to better understand and possibly manage processes 
> being studied.  I will offer some thoughts on this, and I hope if there 
> are errors in my logic or if someone else has a better idea, we will 
> both benefit from their comments.
> 
>   The ideal might be an estimate of something like a mean square for 
> a particular effect to compare with an estimated variance component.
> Such mean squares were a mandatory component of any analysis of variance 
> table prior to the (a) popularization of generalized linear models and 
> (b) availability of software that made it feasible to compute maximum 
> likelihood estimates routinely for unbalanced, mixed-effects models. 
> However, anova(lme(...)) such mean squares are for most purposes 
> unnecessary cluster in a modern anova table.
> 
>   To estimate a mean square for a fixed effect, consider the 
> following log(likelihood) for a mixed-effects model:
> 
>   lglk = (-0.5)*(n*log(2*pi*var.e)-log(det(W)) + 
> t(y-X%*%b)%*%W%*%(y-X%*%b)/var.e),
> 
> where n = the number of observations,
> 
>   b = the fixed-effect parameter variance,
> 
> and the covariance matrix of the residuals, after integrating out the 
> random effects is var.e*solve(W).  In this formulation, the matrix "W" 
> is a function of the variance components.  Since they are not needed to 
> compute the desired mean squares, they are suppressed in the notation here.
> 
>   Then, the maximum likelihood estimate of
> 
>   var.e = SSR/n,
> 
> where SSR = t(y-X%*%b)%*%W%*%(y-X%*%b).
> 
>   Then
> 
>   mle.lglk = (-0.5)*(n*(log(2*pi*SSR/n)-1)-log(det(W))).
> 
>   Now let
> 
>   SSR0 = this generalized sum of squares of residuals (SSR) without 
> effect "1",
> 
> and
> 
>   SSR1 = this generalized SSR with this effect "1".
> 
>   If I've done my math correctly, then
> 
>   D = deviance = 2*log(likelihood ratio)
> = (n*log(SSR0/SSR1)+log(det(W1)/det(W0)))
> 
>   For roughly half a century, a major part of "the analysis of 
> variance" was the Pythagorean idea that the sum of squares under H0 was 
> the sum of squares under H1 plus the sum of squares for effect "1":
> 
>   SSR0 = SS1 + SSR1.
> 
>   Whence,
> 
>   exp((D/n)-log(det(W1)/det(W0))) = 1+SS1/SSR1.
> 
> Thus,
> 
>   SS1 = SSR1*(exp((D/n)-log(det(W1)/det(W0)))-1).
> 
>   If the difference between deg(W1) and det(W0) can be ignored, we get:
> 
>   SS1 = SSR1*(exp((D/n)-1).
> 
>   Now compute MS1 = SS1/df1, and compare with the variance components.
> 
>   If there is a flaw in this logic, I hope someone will disabuse me 
> of it.
> 
>   If this seems too terse or convoluted to follow, please provide a 
> simple, self-contained example, as suggested in the posting guide! 
> "www.R-project.org/posting-guide.html".  You asked a theoretical 
> question, you got a theoretical answer.  If you want a concrete answer, 
> it might help to pose a more concrete question.
> 
>   Hope this helps.
>   Spencer Graves   
> 
> Bruno L. Giordano wrote:
>> Hello,
>> Is there a way to compare the relative relevance of fixed and random 
>> effects in mixed models? I have in mind measures of effect size in 
>> ANOVAs, and would like to obtain similar information with mixed models.
>>
>> Are there information criteria that allow to compare the relevance of 
>> each of the effects in a mixed model to the overall fit?
>>
>> 

Re: [R] Effect size in mixed models

[Bruno L. Giordano]

OK and thanks!! I will try to apply the approach you suggest to my data.

I found a possibly interesting reference on this topic (if my fast reading 
is error-free):

E.  Demidenko, 2004, Mixed Models: Theory and Applications (John Wiley & 
Sons), Chapter 9 (Diagnostics and influence analysis)

outlines a "perturbation-based" methodology called Infinitesimal influence 
analysis

which analyzes the stability of the solution as a function of perturbations 
in the predictors (or in the isolated observations). In this case, 
predictors whose perturbation leads to greater changes in the solution are 
assumed to be more relevant. In the same chapter the idea is extended to 
mixed models, although it focuses on the influence of clusters of 
observations on the solution, and not on the influence of random effects.

A really rough idea I had, more in line with my current mathematical 
background, is to simply examine the percent increase in the residual 
variance when one (or more) fixed/random effect(s) is(are) removed from the 
model.

With this brute approach I found, for example, that for the outcome I am 
analyzing (4 random + 4 fixed effects with interactions) the entire fixed 
portion is less relevant than any of the random effects in isolation.

Bruno



- Original Message - 
From: "Spencer Graves" <[EMAIL PROTECTED]>
To: "Bruno L. Giordano" <[EMAIL PROTECTED]>
Cc: 
Sent: Monday, June 19, 2006 10:52 PM
Subject: Re: [R] Effect size in mixed models


>   You have asked a great question:  It would indeed be useful to compare 
> the relative magnitude of fixed and random effects, e.g. to prioritize 
> efforts to better understand and possibly manage processes being studied. 
> I will offer some thoughts on this, and I hope if there are errors in my 
> logic or if someone else has a better idea, we will both benefit from 
> their comments.
>
>   The ideal might be an estimate of something like a mean square for a 
> particular effect to compare with an estimated variance component.
> Such mean squares were a mandatory component of any analysis of variance 
> table prior to the (a) popularization of generalized linear models and (b) 
> availability of software that made it feasible to compute maximum 
> likelihood estimates routinely for unbalanced, mixed-effects models. 
> However, anova(lme(...)) such mean squares are for most purposes 
> unnecessary cluster in a modern anova table.
>
>   To estimate a mean square for a fixed effect, consider the following 
> log(likelihood) for a mixed-effects model:
>
>   lglk = (-0.5)*(n*log(2*pi*var.e)-log(det(W)) + 
> t(y-X%*%b)%*%W%*%(y-X%*%b)/var.e),
>
> where n = the number of observations,
>
>   b = the fixed-effect parameter variance,
>
> and the covariance matrix of the residuals, after integrating out the 
> random effects is var.e*solve(W).  In this formulation, the matrix "W" is 
> a function of the variance components.  Since they are not needed to 
> compute the desired mean squares, they are suppressed in the notation 
> here.
>
>   Then, the maximum likelihood estimate of
>
>   var.e = SSR/n,
>
> where SSR = t(y-X%*%b)%*%W%*%(y-X%*%b).
>
>   Then
>
>   mle.lglk = (-0.5)*(n*(log(2*pi*SSR/n)-1)-log(det(W))).
>
>   Now let
>
>   SSR0 = this generalized sum of squares of residuals (SSR) without effect 
> "1",
>
> and
>
>   SSR1 = this generalized SSR with this effect "1".
>
>   If I've done my math correctly, then
>
>   D = deviance = 2*log(likelihood ratio)
> = (n*log(SSR0/SSR1)+log(det(W1)/det(W0)))
>
>   For roughly half a century, a major part of "the analysis of variance" 
> was the Pythagorean idea that the sum of squares under H0 was the sum of 
> squares under H1 plus the sum of squares for effect "1":
>
>   SSR0 = SS1 + SSR1.
>
>   Whence,
>
>   exp((D/n)-log(det(W1)/det(W0))) = 1+SS1/SSR1.
>
> Thus,
>
>   SS1 = SSR1*(exp((D/n)-log(det(W1)/det(W0)))-1).
>
>   If the difference between deg(W1) and det(W0) can be ignored, we get:
>
>   SS1 = SSR1*(exp((D/n)-1).
>
>   Now compute MS1 = SS1/df1, and compare with the variance components.
>
>   If there is a flaw in this logic, I hope someone will disabuse me of it.
>
>   If this seems too terse or convoluted to follow, please provide a 
> simple, self-contained example, as suggested in the posting guide! 
> "www.R-project.org/posting-guide.html".  You asked a theoretical question, 
> you got a theoretical answer.  If you want a concrete answer, it might 
> help to pose a more concrete question.
>
>   Hope this helps.
>   Spencer Graves
> Bruno L. Giordano wrote:
>> Hello,
>>

Re: [R] Effect size in mixed models

[Spencer Graves]

  You have asked a great question:  It would indeed be useful to 
compare the relative magnitude of fixed and random effects, e.g. to 
prioritize efforts to better understand and possibly manage processes 
being studied.  I will offer some thoughts on this, and I hope if there 
are errors in my logic or if someone else has a better idea, we will 
both benefit from their comments.

  The ideal might be an estimate of something like a mean square for a 
particular effect to compare with an estimated variance component.
Such mean squares were a mandatory component of any analysis of variance 
table prior to the (a) popularization of generalized linear models and 
(b) availability of software that made it feasible to compute maximum 
likelihood estimates routinely for unbalanced, mixed-effects models. 
However, anova(lme(...)) such mean squares are for most purposes 
unnecessary cluster in a modern anova table.

  To estimate a mean square for a fixed effect, consider the following 
log(likelihood) for a mixed-effects model:

  lglk = (-0.5)*(n*log(2*pi*var.e)-log(det(W)) + 
t(y-X%*%b)%*%W%*%(y-X%*%b)/var.e),

where n = the number of observations,

   b = the fixed-effect parameter variance,

and the covariance matrix of the residuals, after integrating out the 
random effects is var.e*solve(W).  In this formulation, the matrix "W" 
is a function of the variance components.  Since they are not needed to 
compute the desired mean squares, they are suppressed in the notation 
here.

  Then, the maximum likelihood estimate of

  var.e = SSR/n,

where SSR = t(y-X%*%b)%*%W%*%(y-X%*%b).

  Then

  mle.lglk = (-0.5)*(n*(log(2*pi*SSR/n)-1)-log(det(W))).

  Now let

  SSR0 = this generalized sum of squares of residuals (SSR) without 
effect "1",

and

  SSR1 = this generalized SSR with this effect "1".

  If I've done my math correctly, then

  D = deviance = 2*log(likelihood ratio)
= (n*log(SSR0/SSR1)+log(det(W1)/det(W0)))

  For roughly half a century, a major part of "the analysis of 
variance" was the Pythagorean idea that the sum of squares under H0 was 
the sum of squares under H1 plus the sum of squares for effect "1":

  SSR0 = SS1 + SSR1.

  Whence,

  exp((D/n)-log(det(W1)/det(W0))) = 1+SS1/SSR1.

Thus,

  SS1 = SSR1*(exp((D/n)-log(det(W1)/det(W0)))-1).

  If the difference between deg(W1) and det(W0) can be ignored, we get:

  SS1 = SSR1*(exp((D/n)-1).

  Now compute MS1 = SS1/df1, and compare with the variance components.

  If there is a flaw in this logic, I hope someone will disabuse me of 
it.

  If this seems too terse or convoluted to follow, please provide a 
simple, self-contained example, as suggested in the posting guide! 
"www.R-project.org/posting-guide.html".  You asked a theoretical 
question, you got a theoretical answer.  If you want a concrete answer, 
it might help to pose a more concrete question.

  Hope this helps.
  Spencer Graves

Bruno L. Giordano wrote:
> Hello,
> Is there a way to compare the relative relevance of fixed and random effects 
> in mixed models? I have in mind measures of effect size in ANOVAs, and would 
> like to obtain similar information with mixed models.
> 
> Are there information criteria that allow to compare the relevance of each 
> of the effects in a mixed model to the overall fit?
> 
> Thank you,
> Bruno
> 
> __
> R-help@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@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