[R] Estimates with lme(...varPower())
Douglas Bates
bates at stat.wisc.edu
Tue Apr 29 16:32:49 CEST 2003
"Dieter Menne" <dieter.menne at menne-biomed.de> writes:
> we have a 2*3 factorial design, where 2 out of the 3 treatment levels
> are tested on each subject on two different days, each before and after a
> meal (When)
>
> There is strong evidence for heteroscedascticity. The lme-analysis with
> varPower-weighting has significantly lower AIC, estimated power is 1.33.
>
> Estimated values are physiologically reasonable and close to averages for
> the unweighted case (e.g. 130.3+184.8=315 for WhenPost at the base level),
> but for the weighted case they are not (110.2+134.9=245).
>
> Does using weighted estimates make sense at all? Or did I misuse weighted
> lme here? Is there any warning signal I overlooked?
>
> Dieter Menne
>
> -------------------
> No weigthing
>
> > vTonus.lme<-lme(vTonus~Treat*When,data=to,random=~1|Subj/Day)
> AIC BIC logLik
> 844 864 -413
> ....
> Value Std.Error DF t-value p-value
> (Intercept) 130.3 35.3 33 3.69 0.0008
> TreatCel 62.3 45.4 16 1.37 0.1890
> TreatDic 63.4 45.4 16 1.40 0.1818
> WhenPost 184.8 34.8 33 5.32 <.0001
> TreatCel:WhenPost -59.8 49.1 33 -1.22 0.2320
> TreatDic:WhenPost -68.4 49.1 33 -1.39 0.1735
> Correlation:
> (Intr) TretCl TretDc WhnPst TrC:WP
> TreatCel -0.642
> TreatDic -0.642 0.500
> WhenPost -0.492 0.383 0.383
> TreatCel:WhenPost 0.348 -0.541 -0.271 -0.707
> TreatDic:WhenPost 0.348 -0.271 -0.541 -0.707 0.500
>
> With vaPower weighting:
>
> >vTonusW.lme<-lme(vTonus~Treat*When,data=to,
> + weights=varPower(form=~vTonus),random=~1|Subj/Day)
>
> AIC BIC logLik
> 830 852 -405
> ..
> Parameter estimates:
> power 1.33
> ..
> Value Std.Error DF t-value p-value
> (Intercept) 110.2 20.9 33 5.27 <.0001
> TreatCel 54.0 28.5 16 1.90 0.0762
> TreatDic 31.9 28.2 16 1.13 0.2744
> WhenPost 134.9 26.2 33 5.15 <.0001
> TreatCel:WhenPost -84.4 28.7 33 -2.94 0.0060
> TreatDic:WhenPost -96.5 32.9 33 -2.93 0.0061
> Correlation:
> (Intr) TretCl TretDc WhnPst TrC:WP
> TreatCel -0.562
> TreatDic -0.554 0.386
> WhenPost -0.125 0.093 0.091
> TreatCel:WhenPost 0.111 -0.119 -0.078 -0.912
> TreatDic:WhenPost 0.096 -0.070 -0.195 -0.794 0.724
It would be more common to use the predicted response rather than the
observed response for the varPower weights, as in
varPower(form = ~fitted(.))
If you use the observed responses then abnormally low observed
responses will have unusually high weights in the fit, leading to low
predictions, as you observed.
More information about the R-help
mailing list