[R] [R ] Query : problems with the arithmetic operator "^" with function "lme"
MARTIN Ludovic
martinl at mathinfo.ens.univ-reims.fr
Thu May 22 16:25:58 CEST 2003
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.
More information about the R-help
mailing list