1. This is not an R question; it is a statistical issue.
2. R-sig-mixed-models is the appropriate list, not r-help.
-- Bert
On Tue, Jun 26, 2012 at 3:28 AM, KL wrote:
> Hi,
>
> I performed an experiment where I raised different families coming from two
> different source populations, where each family was split up into a
> different treatments. After the experiment I measured several traits on
> each
> individual.
> To test for an effect of either treatment or source as well as their
> interaction, I used a linear mixed effect model with family as random
> factor, i.e.
> lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")
>
> so far so good,
> Now I have to calculate the relative variance components, i.e. the
> percentage of variation that is explained by either treatment or source as
> well as the interaction.
>
> Without a random effect, I could easily use the sums of squares (SS) to
> calculate the variance explained by each factor. But for a mixed model
> (with
> ML estimation), there are no SS, hence I thought I could use Treatment and
> Source as random effects too to estimate the variance, i.e.
>
> lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")
>
> However, in some cases, lme does not converge, hence I used lmer from the
> lme4 package:
>
> lmer(Trait~1+(Treatment*Source|Family),data=DATA)
>
> Where I extract the variances from the model using the summary function:
>
> model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
> results<-model@REmat
> variances<-results[,3]
>
> I get the same values as with the VarCorr function. I use then these values
> to calculate the actual percentage of variation taking the sum as the total
> variation.
>
> Where I am struggling is with the interpretation of the results from the
> initial lme model (with treatment and source as fixed effects) and the
> random model to estimate the variance components (with treatment and source
> as random effect). I find in most cases that the percentage of variance
> explained by each factor does not correspond to the significance of the
> fixed effect.
>
> For example for the trait HD,
> The initial lme suggests a tendency for the interaction as well as a
> significance for Treatment. Using a backward procedure, I find that
> Treatment has a close to significant tendency. However, estimating variance
> components, I find that Source has the highest variance, making up to 26.7%
> of the total variance.
>
>
>
> anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
> numDF denDF F-value p-value
>(Intercept)1 426 0.044523 0.8330
>as.factor(Treatment) 1 426 5.935189 0.0153
>as.factor(Source) 111 0.042662 0.8401
>as.factor(Treatment):as.factor(Source) 1 426 3.754112 0.0533
>
>
>
>
>
>
>
> summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
>Linear mixed model fit by REML
>Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family)
> Data: regrexpdat
>AICBIC logLik deviance REMLdev
> -103.5 -54.43 63.75 -132.5 -127.5
>Random effects:
> Groups Name Variance Std.Dev.
> Corr
> Family (Intercept) 0.0113276 0.106431
> as.factor(Treatment) 0.0063710 0.079819
> 0.405
> as.factor(Source) 0.0235294 0.153393
> -0.134 -0.157
> as.factor(Treatment)L:as.factor(Source) 0.0076353 0.087380
> -0.578 -0.589 -0.585
> Residual 0.0394610 0.198648
>Number of obs: 441, groups: Family, 13
>
>Fixed effects:
>Estimate Std. Error t value
>(Intercept) -0.027400.03237 -0.846
>
>
>
> Hence my question is, is it correct what I am doing? Or should I use
> another
> way to estimate the amount of variance explained by each factor (i.e.
> Treatment, Source and their interaction). For example, would the effect
> sizes be a more appropriate way to go?
>
>
> Thanks!
>
> Kay Lucek
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-estimate-variance-components-with-lmer-for-models-with-random-effects-and-compare-them-with-ls-tp4634492.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Bert Gunter
Genentech Nonclinical Biostatistics
Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-
