R help - May 2003 - [R ] Query : problems with the arithmetic operator "^" with function "lme"

Dear all,

I've got a problem in including square variables in lme function. I've 
tried to work on Dialyzer data of Pinheiro and Bates'book. 

We fit the heteroscedastic model with:
> data(Dialyzer)
> fm2Dial.lme<-lme(rate~(pressure+pressure^2+pressure^3+pressure^4)*QB,
+    Dialyzer,~pressure+pressure^2,weights=varPower(form=~pressure))

We Obtain
> fm2Dial.lme
Linear mixed-effects model fit by REML
  Data: Dialyzer 
  Log-restricted-likelihood: -488.4535
  Fixed: rate ~ (pressure + pressure^2 + pressure^3 + pressure^4) * QB 
   (Intercept)       pressure          QB300 pressure:QB300 
     39.362769       1.480331      -4.547449       7.515690 

Random effects:
 Formula: ~pressure + pressure^2 | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 0.01102331 (Intr)
pressure    1.26972080 -0.823
Residual    8.79053329       

Variance function:
 Structure: Power of variance covariate
 Formula: ~pressure 
 Parameter estimates:
    power 
-1.025219 
Number of Observations: 140
Number of Groups: 20 

They are not coefficients associated with the "pressure^2, 
pressure^3, ..."
However, the model called is " rate ~ (pressure + pressure^2 + 
pressure^3 + pressure^4) * QB " . "^" is a problem !
So, we fit the model like this, including two matrices, for the fixed 
effects and the random effects:
>Dialyzer$PressureA<-cbind(Dialyzer$pressure,...,Dialyzer$pressure^4)
>Dialyzer$PressureB<-cbind(Dialyser$pressure,Dialyzer$pressure^2)
Now, we fit the same model like this:
>fm3Dial.lme<-lme(rate~(PressureA)*QB,
+    Dialyzer,~PressureB,weights=varPower(form=~pressure))

We obtain:

Linear mixed-effects model fit by REML
  Data: Dialyzer 
  Log-restricted-likelihood: -309.5058
  Fixed: rate ~ (PressureA) * QB 
   (Intercept)     PressureA1     PressureA2     PressureA3     
   -17.6805986     93.7145527    -49.1906874     12.2471222      
. . . 
Random effects:
 Formula: ~PressureB | Subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr         
(Intercept) 1.859195 (Intr) PrssB1
PressureB1  5.327444 -0.523       
PressureB2  1.648356  0.364 -0.954
Residual    1.261931              
. . .
 
Now, it's OK. The anova method returns
> anova(fm3Dial.lme)
             numDF denDF  F-value p-value
(Intercept)      1   112 551.5492  <.0001
PressureA        4   112 968.6231  <.0001
QB               1    18   4.7268  0.0433
PressureA:QB     4   112  20.9273  <.0001

However, the anova method is used to test the significanse of the terms 
in the order they were entered in the model. In Pinheiro and 
Bates'book, the result is
>anova(fm2Dial.lme)
 
                   numDF denDF  F-value p-value
      (Intercept)      1   112    552.9  <.0001
         pressure      1   112   2328.6  <.0001
    I(pressure^2)      1   112   1174.6  <.0001
           ...         ... ...     ...   ...
  I(pressure^2):QB     1   112      1.4  0.2477
  I(pressure^3):QB     1   112      2.2  0.1370
  I(pressure^4):QB     1   112      0.2  0.6840

The three large p-values suggest they are not needed in the model. They 
could be elimated from the model. It isn't indicated when we use 
matrices PressureA and PressureB ! So, how can we do? 
  
I would be grateful if anyone could help me.

Cordially,

Martin Ludovic.

"MARTIN ?Ludovic" <martinl at mathinfo.ens.univ-reims.fr>
writes:
> > data(Dialyzer)
> >
fm2Dial.lme<-lme(rate~(pressure+pressure^2+pressure^3+pressure^4)*QB,
> +    Dialyzer,~pressure+pressure^2,weights=varPower(form=~pressure))
You'll need I(pressure^2) in R. (FAQ 3.3.2)


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907