[R] contrasts in lm and lme

[Peter B. Mandeville]

I am using RW 1.2.3. on an IBM PC 300GL.

Using the data bp.dat which accompanies

    Helen Brown and Robin Prescott
    1999 Applied Mixed Models in Medicine. Statistics in Practice. 
         John Wiley & Sons, Inc., New York, NY, USA

which is also found at www.med.ed.ac.uk/phs/mixed. The data file was opened
and initialized with

> dat <- read.table("bp.dat")
> names(dat) <-
c("patient","visit","center","treatment","dbp","dbp1","cf","cf1")
> attach(dat)
> Patient <- factor(patient)
> Treatment <- factor(treatment)
> Center <- factor(center)
> Visit <- factor(visit)
> dat1 <- data.frame(Patient,Visit,Center,Treatment,dbp,dbp1)
> sapply(dat1,data.class)
  Patient     Visit    Center Treatment       dbp      dbp1 
 "factor"  "factor"  "factor"  "factor" "numeric" "numeric" 

After which the following code was run

> library(nlme) # needed for the contrast contr.SAS
> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> res <- lm(dbp~Treatment+Visit+dbp1)
> anova(res)
> summary(res)

and was repeated leaving out library(nlme) and replacing

> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))

with the following contrasts

> options(contrasts=c(factor="contr.treatment",ordered="contr.poly"))
> options(contrasts=c(factor="contr.helmert",ordered="contr.poly"))
> options(contrasts=c(factor="contr.sum",ordered="contr.poly"))
> library(MASS) # needed for the contrast contr.sdif
> options(contrasts=c(factor="contr.sdif",ordered="contr.poly"))

The results from anova were igual. For example for the factor Treatment had
the same probabilities in every case but

         SAS        treatment    helmert    sum        sdif
Pr(>F)   3.073e-05  3.073e-05    3.073e-05  3.073e-05  3.073e-05

the probabilities of the different contrasts were

         SAS        treatment    helmert    sum        sdif
A-B                 0.072635     0.072635              0.072635
A-C      8.69e-06   8.69e-06                0.000322
B-C      0.00875                            0.637971   0.00875

Why does the contrast contr.sum have distinct results from the other
contrasts? Which ones are confiable?

Pinheiro and Bates 2000:17 state that

    Although the individual parameter estimates for the Type factor are
different
    between the two fits, the anova resultas are the same. The difference
in the
    parameter estimates simply reflects the fact that different contrasts
are being
    estimated.

If the process is repeated with lme in place of lm with

> res <- lm(dbp~Treatment+Visit+dbp1,random=~1|Patient)

in place of 

> res <- lm(dbp~Treatment+Visit+dbp1)


         SAS        treatment    helmert    sum        sdif
AIC      7501.212   7501.212     7511.151   7506.182   7501.212
BIC      7546.116   7546.116     7556.055   7551.086   7546.116
logLik   -3741.606  -3741.606    -3746.576  -3744.091  -3741.606

Given that AIC and BIC are calculated logLik, it is reasonable that they
differ given the different values of the logLik, but is it reasonable that
the logLik's are different?

         SAS        treatment    helmert    sum        sdif
A-B                 0.2330       0.2330                0.2330
A-C      0.0033     0.0033                  0.0168
B-C      0.0850                             0.7553     0.0850

Again, why does the contrast contr.sum have distinct results from the other
contrasts? Which ones are confiable?

Thank you very much

Peter B.


--
Peter B. Mandeville                             mandevip at deimos.tc.uaslp.mx
Jefe del Depto. de Informática y Bioestadística rpe1531 at pasteur.fmed.uaslp.mx 
Facultad de Medicine                            Tel: 48 26-23-45 ext. 232
Universidad Autónoma de San Luis Potosí         Fax: 48 28-23-52
Av. V. Carranza 2405
Col. Los Filtros
Apartado Postal 145
San Luis Potosí, S.L.P.
78210 México

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._