Mixed-effects heterocedastic model

[artbio]

I fitted a nixed-effects model to a dataset from a longitudinal study in R. This is the R code I used:

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:

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.

 

[mmwave]

Thank you for sharing your code and output with us. It seems like you have fit a mixed-effects model with a heteroscedastic structure, meaning that the variance of your outcome variable is different for each level of the categorical variable "lv". This is a common approach in longitudinal studies, where the variability of the outcome may change over time or across different groups.

To express this model in mathematical terms, you can use the following equation:

Y_ij = β_0 + β_1*X_1 + β_2*X_2 + ... + β_p*X_p + b_0i + b_1i*X_1 + b_2i*X_2 + ... + b_qi*X_q + e_ij

Where:
- Y_ij is the observed outcome for individual i at time j
- β_0 is the overall intercept
- β_1 to β_p are the fixed effects coefficients for the independent variables X_1 to X_p
- b_0i is the random intercept for individual i
- b_1i to b_qi are the random slope coefficients for individual i
- e_ij is the residual error term for individual i at time j

The key difference between this equation and a regular mixed-effects model is the inclusion of b_0i and b_1i, which represent the random effects for the intercept and slope, respectively. These random effects are modeled as having a normal distribution with mean 0 and variance σ^2_b0 and σ^2_b1, respectively. In your case, you have specified that the variance of these random effects varies by the level of the categorical variable "lv", which is why you have included the weights=varIdent(form=~1|lv) argument in your code.

I hope this helps clarify the mathematical representation of your heteroscedastic mixed-effects model. If you have any further questions, please don't hesitate to ask.