纵向数据

I've been working with the R Orthodont dataset in the "nlme" package. Just use install.packages("nlme");library(nlme);head(Orthodont) to take a look. The dataset is comprised of distance between the pituitary and the pterygomaxillary fissure measured in 27 children over time. Using the lme4 package I can fit a nonlinear mixed effects model using a logistic curve as my functional form. I can choose to have the asymptote and midpoint entered as random effects

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 包是否可以做到这一点？

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

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

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

# 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))

发生的事情是当性别==0（男），模型实现数值：

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

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

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

因此，随着年龄的增长，我们达到的渐近水平实际上是不对称+不对称2， 不只是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。如果您想提取非线性固定效应，则模型的可视化需要进行一些修改（与中线性模型的可视化类似）如何通过观察提取 lmer 固定效应？ ):

# 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似乎有点棘手。