[R] Estimates with lme(...varPower())

[Douglas Bates]

"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.