[R-sig-ME] Poor fits with corCAR1 in lme and phi estimate near 1

[Weigand, Stephen D.]

I'm trying to understand why I get poor fits with a lme model
using CAR1 autocorrelation structure.

I'm modeling the volume of a part of a person's brain over time 
using lme. The volume depends on head size, sex, age, and age^2.
I give each person a random intercept.

To give you an idea of my sample sizes here is

table(table(mydata$group))

  3  4  5  6  7  8
16 16  4  7  1  2

My lme call is like this

   ccc <- list(returnObject = TRUE,
               msVerbose = TRUE,
               maxIter = 200,
               msMaxIter = 200,
               msMaxEval = 500)

   fit <-
     lme(vol ~ size + sex + age + I(age^2),
         random = ~ 1 | id,
         correlation = corCAR1(form = ~ age | id),
         data = dG1,
         control = ccc)

I get convergence in 49 iterations with no warnings.

My first indication of a problem is that I get an estimate of phi 
to be 0.996. What's happening is that the fitted values are all 
above or below the observed values. Here's the data for four 
people with a call to xyplot() below to illustrate.

mydata <- structure( list( y = c(101.513696992828,
108.67123416, 121.77207312, 121.364514, 122.44394352,
126.28247856, 66.4450535719635, 71.06907726, 73.34944657,
77.73755495, 73.712, 77.188, 77.802, 22.6576319594522,
24.08128381, 23.3368159, 25.57910501, 26.39909865,
27.52373389 ), age = c(81.4784394250513, 83.4469541409993,
86.7405886379192, 89.0184804928131, 90.3655030800821,
92.413415468857, 73.1663244353183, 76.0602327173169,
78.1629021218344, 79.7125256673511, 80.807665982204,
82.9075975359343, 84.1697467488022, 84.3997262149213,
85.719370294319, 87.0581793292266, 88.558521560575,
91.419575633128, 93.5277207392197 ), g = structure(c(4, 4,
4, 4, 4, 4, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3), .Label =
c("1", "2", "3", "4"), class = "factor"), yhat =
c(58.5700637637191, 61.2504094310121, 65.7874577250516,
68.9636770643335, 70.8566873027208, 73.7557021229065,
55.8860711627563, 59.8790531221204, 62.8120620013217,
64.9907426560824, 71.0988108284064, 74.0280180708455,
75.8014293467406, 50.2485578022043, 52.0817661160286,
53.9523602784579, 56.0615306867955, 60.1213127021499,
63.1444086154206 )), .Names = c("y", "age", "g", "yhat"),
row.names = c(NA, 19), class = "data.frame")

### plot illustrating the "missed fits":
require(lattice)
xyplot(y + yhat ~ age | g, data = mydata,
        xlab = "Age",
        ylab = "Volume",
        lwd = 2,
        type = "b",
        col.symbol = c("black", "transparent"),
        col.line = c("transparent", "black"))


When I use Gaussian-decay autocorrelation, things look fine. I'm 
wondering if anybody has insight into this. I'd be happy to stay 
away from the CAR1 autocorrelation but in similar models the 
likelihood is higher with CAR1 than with Gaussian decay (even 
when there are no problems with the CAR1 model).

Thanks for any thoughts,

Stephen


-- 
::::::::::::::::::::::::::::::::::
Stephen Weigand
Division of Biostatistics
Mayo Clinic Rochester, Minn., USA
Phone (507) 266-1650, fax 284-9542