纵向数据

Question

我一直在使用“nlme”包中的 R Orthodont 数据集。只需使用install.packages("nlme");library(nlme);head(Orthodont) 即可查看。该数据集包括随时间推移在 27 名儿童中测量的垂体和翼上颌裂之间的距离。 使用 lme4 包，我可以使用逻辑曲线作为我的函数形式来拟合非线性混合效应模型。我可以选择将渐近线和中点作为随机效应输入

nm1 <- nlmer(distance ~ SSlogis(age,Asym, xmid, scal) ~ (Asym | Subject) + (xmid | Subject), Orthodont, start = c(Asym =25,xmid = 11, scal = 3), corr = FALSE,verb=1)

我真正想知道的是性别是否会改变这些参数。不幸的是，在线示例不包括主题和组示例。这甚至可以使用 lme4 包吗？

Answer 1

我相信可以通过为自定义模型公式及其渐变创建函数来做到这一点。标准 SSlogis 函数使用以下形式的逻辑函数：

f(input) = Asym/(1+exp((xmid-input)/scal)) # as in ?SSlogis

您可以修改上述语句以满足您的需要，而不是调用 SSlogis。我相信你会想看看性别是否对固定效应有影响。这是在 Asym2 中修改特定性别的 Asym 子群体效果的示例代码：

# Just for loading the data, we will use lme4 for model fitting, not nlme
library(nlme)
library(lme4)
# Careful when loading both nlme and lme4 as they have overlap, strange behaviour may occur

# A more generalized form could be taken e.g. from http://en.wikipedia.org/wiki/Generalised_logistic_curve
# A custom model structure:
Model <- function(age, Asym, Asym2, xmid, scal, Gender) 
{
    # Taken from ?SSlogis, standard form:
    #Asym/(1+exp((xmid-input)/scal))
    # Add gender-specific term to Asym2
    (Asym+Asym2*Gender)/(1+exp((xmid-age)/scal))
    # Evaluation of above form is returned by this function
}

# Model gradient, notice that we include all 
# estimated fixed effects like 'Asym', 'Asym2', 'xmid' and 'scal' here,
# but not covariates from the data: 'age' and 'Gender'
ModelGradient <- deriv(
    body(Model)[[2]], 
    namevec = c("Asym", "Asym2", "xmid", "scal"), 
    function.arg=Model
)

引入性别效应的一种相当典型的方法是使用二进制编码。我会将 Sex 变量转换为二进制编码的Gender：

# Binary coding for the gender
Orthodont2 <- data.frame(Orthodont, Gender = as.numeric(Orthodont[,"Sex"])-1)
#> table(Orthodont2[,"Gender"])
# 0  1 
#64 44 
# Ordering data based on factor levels so they don't mix up paneling in lattice later on
Orthodont2 <- Orthodont2[order(Orthodont2[,"Subject"]),]

然后我可以拟合自定义模型：

# Fit the non-linear mixed effects model
fit <- nlmer(
    # Response
    distance ~ 
    # Fixed effects
    ModelGradient(age = age, Asym, Asym2, xmid, scal, Gender = Gender) ~ 
    # replaces: SSlogis(age,Asym, xmid, scal) ~ 
    # Random effects
    (Asym | Subject) + (xmid | Subject), 
    # Data
    data = Orthodont2, 
    start = c(Asym = 25, Asym2 = 15, xmid = 11, scal = 3))

当Gender==0（男性）时，模型实现了价值：

(Asym+Asym2*0)/(1+exp((xmid-age)/scal)) = (Asym)/(1+exp((xmid-age)/scal))

这实际上是标准的 SSlogis 函数形式。但是，现在有一个二进制开关，如果 Gender==1（女性）：

(Asym+Asym2)/(1+exp((xmid-age)/scal))

所以我们随着年龄的增长达到的渐近水平实际上是 Asym + Asym2，而不仅仅是 Asym，对于女性个体。

还请注意，我没有为 Asym2 指定新的随机效果。由于 Asym 对性别没有特异性，因此女性个体也可能因 Asym 项而在其个体渐近水平上存在差异。模型拟合：

> summary(fit)
Nonlinear mixed model fit by the Laplace approximation 
Formula: distance ~ ModelGradient(age = age, Asym, Asym2, xmid, scal,      Gender = Gender) ~ (Asym | Subject) + (xmid | Subject) 
   Data: Orthodont2 
   AIC   BIC logLik deviance
 268.7 287.5 -127.4    254.7
Random effects:
 Groups   Name Variance Std.Dev.
 Subject  Asym 7.0499   2.6552  
 Subject  xmid 4.4285   2.1044  
 Residual      1.5354   1.2391  
Number of obs: 108, groups: Subject, 27

Fixed effects:
      Estimate Std. Error t value
Asym    29.882      1.947  15.350
Asym2   -3.493      1.222  -2.859
xmid     1.240      1.068   1.161
scal     5.532      1.782   3.104

Correlation of Fixed Effects:
      Asym   Asym2  xmid  
Asym2 -0.471              
xmid  -0.584  0.167       
scal   0.901 -0.239 -0.773

看起来可能存在性别特定效应 (t -2.859)，因此随着“年龄”的增加，女性患者的“距离”值似乎会降低：29.882 - 3.493 = 26.389

我并不一定建议这是一个好的/最佳模型，只是展示了您可以如何继续自定义 lme4 中的非线性模型。如果您想提取非线性固定效应，模型的可视化需要一些修改（与How do I extract lmer fixed effects by observation? 中线性模型的可视化类似）：

# Extracting fixed effects components by calling the model function, a bit messy but it works
# I like to do this for visualizing the model fit
fixefmat <- matrix(rep(fixef(fit), times=dim(Orthodont2)[1]), ncol=length(fixef(fit)), byrow=TRUE)
colnames(fixefmat) <- names(fixef(fit))
Orthtemp <- data.frame(fixefmat, Orthodont2)
attach(Orthtemp)
# see str(Orthtemp)
# Evaluate the function for rows of the attached data.frame to extract fixed effects corresponding to observations
fix = as.vector(as.formula(body(Model)[[2]]))
detach(Orthtemp)

nobs <- 4 # 4 observations per subject
legend = list(text=list(c("y", "Xb + Zu", "Xb")), lines = list(col=c("blue", "red", "black"), pch=c(1,1,1), lwd=c(1,1,1), type=c("b","b","b")))
require(lattice)
xyplot(
    distance ~ age | Subject, 
    data = Orthodont2,
    panel = function(x, y, ...){
        panel.points(x, y, type='b', col='blue')
        panel.points(x, fix[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='black')
        panel.points(x, fitted(fit)[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='red')
    },
    key = legend
)

# Residuals
plot(Orthodont2[,"distance"], resid(fit), xlab="y", ylab="e")

# Distribution of random effects
par(mfrow=c(1,2))
hist(ranef(fit)[[1]][,1], xlab="Random 'Asym'", main="")
hist(ranef(fit)[[1]][,2], xlab="Random 'xmid'", main="")
# Random 'xmid' seems a bit skewed to the right and may violate normal distribution assumption
# This is due to M13 having a bit abnormal growth curve (random effects):
#           Asym       xmid
#M13  3.07301310  3.9077583

图形输出：

请注意，上图中女性 (F##) 个体略低于男性 (M##) 个体（黑线）。例如。 M10 F10 中间区域面板的差异。

用于观察指定模型的某些特征的残差和随机效应。单个 M13 似乎有点棘手。