9

我一直在使用“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 包吗？

4

1 回答 1

20

我相信可以通过为自定义模型公式及其梯度创建一个函数来做这样的事情。标准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中的非线性模型。如果要提取非线性固定效应，模型的可视化需要一些修改（与如何通过观察提取 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 似乎有点棘手。

于 2013-06-01T02:18:37.237 回答