R help - Apr 2003 - Estimates with lme(...varPower())

Dear R-List,

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

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

Dear Douglas Bates,
> It would be more common to use the predicted response rather than the
> observed response for the varPower weights, as in
>  varPower(form = ~fitted(.))
Following the example on page 223 of PB, I had used that first, but did not
achieve convergence. Will try with more iterations.
> 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.

Dieter Menne