- #1
artbio
- 14
- 0
I fitted a nixed-effects model to a dataset from a longitudinal study in R. This is the R code I used:
But then I determined that an heterocedastic model with different variances per "lv" level (lv is a categoric variable) would fit my data better. And I fitted this model:
This is the output from R's summary:
My problem is that I don't know exactly how to express this second model in mathematical terms. I am following "Mixed-Effects Models in S and S-Plus" by Pinheiro and Bates. The book seems vague in this regard. Your help would be appreciated.
Thanks.
Code:
data <- read.table("heart1.txt", header=TRUE)
library(MASS)
scope <- list(upper=~sex+age+time+ef+bsa+lvh+prenyha+redo+size+con.cabg+creat+dm+acei+lv+emergenc+hc+sten.reg.mix+hs, lower=~time)
scope$full <- lvmi~sex+age+time+ef+bsa+lvh+prenyha+redo+size+con.cabg+creat+dm+acei+lv+emergenc+hc+sten.reg.mix+hs
scope$reduced <- lvmi~time
library(nlme)
lme5 <- lme(fixed=scope$full, data=data, random=~1+time|num, method="ML")
summary(lme5)
lme6 <- stepAIC(object=lme5, scope=scope, direction="backward")
summary(lme6)
anova(lme5, lme6)
But then I determined that an heterocedastic model with different variances per "lv" level (lv is a categoric variable) would fit my data better. And I fitted this model:
Code:
lme_d <- update(lme, weights=varIdent(form=~1|lv), method="ML")
summary(lme_d)
This is the output from R's summary:
Linear mixed-effects model fit by maximum likelihood
Data: data
AIC BIC logLik
81.31992 142.2811 -26.65996
Random effects:
Formula: ~1 + time | num
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.29954440 (Intr)
time 0.03710599 -0.341
Residual 0.16750748
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | lv
Parameter estimates:
Good Moderate Bad
1.000000 1.187027 1.402825
Fixed effects: lvmi ~ sex + time + bsa + size + lv + hs
Value Std.Error DF t-value p-value
(Intercept) 4.498641 0.3869856 424 11.624827 0.0000
sexFemale -0.165977 0.0657987 143 -2.522504 0.0127
time -0.006066 0.0058608 424 -1.034995 0.3013
bsa -0.357589 0.1155192 143 -3.095495 0.0024
size 0.051367 0.0143853 143 3.570799 0.0005
lvModerate 0.019042 0.0567352 143 0.335625 0.7376
lvBad 0.230576 0.0907100 143 2.541902 0.0121
hsStentless -0.156231 0.0686556 143 -2.275579 0.0244
Correlation:
(Intr) sexFml time bsa size lvMdrt lvBad
sexFemale -0.521
time -0.040 0.027
bsa -0.487 0.349 0.039
size -0.798 0.321 -0.023 -0.123
lvModerate 0.123 -0.079 -0.008 0.003 -0.199
lvBad -0.080 0.104 0.009 -0.016 0.065 0.205
hsStentless 0.384 -0.201 0.017 0.231 -0.650 0.113 -0.087
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.923802823 -0.539808232 0.000763499 0.531272258 3.095870906
Number of Observations: 575
Number of Groups: 150
My problem is that I don't know exactly how to express this second model in mathematical terms. I am following "Mixed-Effects Models in S and S-Plus" by Pinheiro and Bates. The book seems vague in this regard. Your help would be appreciated.
Thanks.