datkins@u.washington.edu
2004-Dec-29 21:55 UTC
[R] Discrepancy between intervals.lme and coef.lme
I'm using R on Windows v2.0.1 with the nlme package (v3.1-53) and am finding
some unexpected discrepancies in the output of intervals.lme and coef.lme.
I've included a toy dataset at the end, but briefly, the data are
longitudinal data from couples in marital therapy. Each spouse's
relationship satisfaction is measured 4 times; I've fit both linear and
quadratic models to the change over time. The quadratic fits show the
discrpancies.
Here's the call to lmList, coef, and intervals:
tmp.lis1 <- lmList(dv ~ time + I(time^2)| id/sex, data = tmp.df,
na.action = na.omit)
coef(tmp.lis1)
intervals(tmp.lis1)
Here is the coef() output:
(Intercept) time I(time^2)
1/Husband 89.60 11.100000 -2.500000
1/Wife 69.80 5.300000 -0.500000
2/Husband 49.00 8.833333 -2.833333
2/Wife 45.00 28.666667 -9.666667
3/Husband 96.00 -6.000000 NA
3/Wife 60.00 19.000000 NA
4/Husband 70.00 48.500000 -16.500000
4/Wife 92.00 14.500000 -4.500000
5/Husband 75.00 43.500000 -14.500000
5/Wife 87.00 37.000000 -14.000000
6/Husband 66.75 1.250000 -1.250000
6/Wife 66.15 12.150000 -2.750000
7/Husband 92.75 6.750000 -0.750000
7/Wife 82.35 17.850000 -3.250000
8/Husband 76.15 -25.350000 11.750000
8/Wife 100.50 -12.000000 6.000000
And just the (Intercept) portion of the intervals() output:
, , (Intercept)
lower est. upper
1/Husband 72.4335719 89.6000000 106.7664281
1/Wife 52.6335719 69.8000000 86.9664281
2/Husband NaN 49.0000000 NaN
2/Wife NaN 45.0000000 NaN
3/Husband NaN 96.0000000 NaN
3/Wife NaN NaN NaN
4/Husband NaN NaN NaN
4/Wife NaN NaN NaN
5/Husband NaN NaN NaN
5/Wife NaN NaN NaN
6/Husband -3.8551591 -0.7453560 2.3644471
6/Wife -3.2306740 -1.4468675 0.3369390
7/Husband -2.1707917 -0.1916630 1.7874658
7/Wife -2.4667397 -1.0766253 0.3134891
8/Husband 4.0996388 4.5693563 5.0390738
8/Wife 0.9368527 1.7888544 2.6408560
Notice that the intercept estimates for couples 6-8 are wildly different between
the coef() and intervals() output. Granted, fitting an intercept, slope, and
quadratic to 4 data points doesn't leave much for an error term, but it
seems like the intercept coefficients ought to be the same. If the quadratic is
dropped, there is no discrepancy between coef() and intevals(), so perhaps this
is related to the complexity of the model vs. sparseness of data?
Any insights appreciated (data below).
Dave
--
Dave Atkins, PhD
datkins at u.washington.edu
tmp.df <- data.frame(id = as.factor(rep(1:8, each=8)),
sex = factor(rep(0:1, each=4, length.out=64), 0:1,
c("Husband","Wife")),
time = rep(0:3, length.out=64),
dv = c(92, 91, 109, 98, 70, 74, 79, 81, 49, 55, NA,
50, 45, 64, NA, 44, NA, 90, 84, NA, NA, 79, 98, NA,
70, 102, 101, NA, 92, 102, 103, NA, 75, 104, 104,
NA, 87, 110, 105, NA, 66, 69, 62, 60, 67, 73, 82,
77, 91, 104, 98, 108, 81, 101, 101, 108, 75, 66,
69, 107, 102, 90, 105, 117))
datkins at u.washington.edu writes:> I'm using R on Windows v2.0.1 with the nlme package (v3.1-53) and am finding some unexpected discrepancies in the output of intervals.lme and coef.lme. I've included a toy dataset at the end, but briefly, the data are longitudinal data from couples in marital therapy. Each spouse's relationship satisfaction is measured 4 times; I've fit both linear and quadratic models to the change over time. The quadratic fits show the discrpancies. > > Here's the call to lmList, coef, and intervals: > > tmp.lis1 <- lmList(dv ~ time + I(time^2)| id/sex, data = tmp.df, > na.action = na.omit) > coef(tmp.lis1) > intervals(tmp.lis1) > > Here is the coef() output: > (Intercept) time I(time^2) > 1/Husband 89.60 11.100000 -2.500000 > 1/Wife 69.80 5.300000 -0.500000 > 2/Husband 49.00 8.833333 -2.833333 > 2/Wife 45.00 28.666667 -9.666667 > 3/Husband 96.00 -6.000000 NA > 3/Wife 60.00 19.000000 NA > 4/Husband 70.00 48.500000 -16.500000 > 4/Wife 92.00 14.500000 -4.500000 > 5/Husband 75.00 43.500000 -14.500000 > 5/Wife 87.00 37.000000 -14.000000 > 6/Husband 66.75 1.250000 -1.250000 > 6/Wife 66.15 12.150000 -2.750000 > 7/Husband 92.75 6.750000 -0.750000 > 7/Wife 82.35 17.850000 -3.250000 > 8/Husband 76.15 -25.350000 11.750000 > 8/Wife 100.50 -12.000000 6.000000 > > And just the (Intercept) portion of the intervals() output: > , , (Intercept) > > lower est. upper > 1/Husband 72.4335719 89.6000000 106.7664281 > 1/Wife 52.6335719 69.8000000 86.9664281 > 2/Husband NaN 49.0000000 NaN > 2/Wife NaN 45.0000000 NaN > 3/Husband NaN 96.0000000 NaN > 3/Wife NaN NaN NaN > 4/Husband NaN NaN NaN > 4/Wife NaN NaN NaN > 5/Husband NaN NaN NaN > 5/Wife NaN NaN NaN > 6/Husband -3.8551591 -0.7453560 2.3644471 > 6/Wife -3.2306740 -1.4468675 0.3369390 > 7/Husband -2.1707917 -0.1916630 1.7874658 > 7/Wife -2.4667397 -1.0766253 0.3134891 > 8/Husband 4.0996388 4.5693563 5.0390738 > 8/Wife 0.9368527 1.7888544 2.6408560 > > Notice that the intercept estimates for couples 6-8 are wildly different between the coef() and intervals() output. Granted, fitting an intercept, slope, and quadratic to 4 data points doesn't leave much for an error term, but it seems like the intercept coefficients ought to be the same. If the quadratic is dropped, there is no discrepancy between coef() and intevals(), so perhaps this is related to the complexity of the model vs. sparseness of data? > > Any insights appreciated (data below).Looks like the sort of things that happens if parts of code don't take notice of pivoting caused by singularities in the model. BTW: This wouldn't look like an issue with the lme methods, but with lmList counterparts. -- 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