如何使用 nlmer() 为重复测量数据拟合非线性混合效应模型?

机器算法验证 r 混合模式 重复测量 lme4-nlme
2022-03-25 10:42:29

我正在尝试分析重复测量数据,并且正在努力使其在R. 我的数据基本上如下,我有两个治疗组。每组中的每个受试者每天都接受测试并给出分数(测试中正确的百分比)。数据为长格式:

Time Percent Subject   Group
   1       0    GK11 Ethanol
   2       0    GK11 Ethanol
   3       0    GK11 Ethanol
   4       0    GK11 Ethanol
   5       0    GK11 Ethanol
   6       0    GK11 Ethanol

数据类似于逻辑曲线,受试者在几天内表现非常糟糕,随后迅速改善,随后进入平台期。我想知道治疗是否对测试性能曲线有影响。我的想法nlmer(). 我可以使用以下方法为每个组拟合线:lme4R

print(nm1 <- nlmer(Percent ~ SSlogis(Time,Asym, xmid, scal) ~ Asym | Subject,
salinedata, start = c(Asym =.60,  xmid = 23, scal = 5)), corr = FALSE)

我可以通过查看估计线的不同参数和标准偏差的估计值来比较组,但我不确定这是不是正确的方法。任何帮助将不胜感激。

1个回答

您可以使用正态似然比检验。这是一个简单的例子。首先,让我们根据您的参数创建 10 个人的观察结果:

Asym = .6
xmid = 23
scal = 5

n = 10
time = seq(1,60,5)

d = data.frame(time=rep(time,10),
               Asym, xmid, scal, group=0)
d$subj = factor(rep(1:n, each=length(time)))

现在让其中一半有不同的渐近线和中点参数:

ind = (nrow(d)/2):nrow(d)
d$Asym[ind] = d$Asym[ind] + .1
d$xmid[ind] = d$xmid[ind] + 10
d$group[ind] = 1
d$group=factor(d$group)

我们可以根据模型模拟所有个体的响应值:

set.seed(1)
d = transform(d, y = Asym/(1+exp((xmid-time)/scal)) +
                     rnorm(nrow(d), sd=.04))
library(lattice)
xyplot(y~time | group, group=subj,
       data=d, type=c("g","l"), col="black")

数据的意大利面条图

我们可以看到两组之间的明显差异,模型应该能够识别的差异。现在让我们首先尝试拟合一个简单的模型,忽略组:

> fm1 = nls(y ~ SSlogis(time, Asym, xmid, scal), data=d)
> coef(fm1)
      Asym       xmid       scal 
 0.6633042 28.5219166  5.8286082

也许正如预期的那样, 和 的估计Asymxmid介于两组的实际参数值之间。(不过,这种情况并不明显,因为比例参数也发生了变化,以调整模型的错误规格。)现在让我们拟合一个完整的模型,两组参数不同:

> fm2 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal[group]),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=rep(5,2)))
> coef(fm2)
    Asym1     Asym2     xmid1     xmid2     scal1     scal2 
 0.602768  0.714199 22.769315 33.331976  4.629332  4.749555

由于这两个模型是嵌套的,我们可以做一个似然比检验:

> anova(fm1, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym, xmid, scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df  Sum Sq F value    Pr(>F)    
1    117    0.70968                                 
2    114    0.13934  3 0.57034  155.54 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

极小的p值清楚地表明简单模型过于简单;两组的参数确实不同。

然而,这两个尺度参数估计值几乎相同,仅相差 0.1。也许我们只需要一个比例参数?(当然我们知道答案是肯定的,因为我们有模拟数据。)

(两个渐近线参数之间的差异也只有 0.1,但当我们考虑标准误差时,差异很大summary(fm2)- 请参阅。)

所以我们拟合了一个新模型,scale两个组有一个共同的参数,但是Asymxmid之前的参数不同:

> fm3 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=5))
> coef(fm3)
     Asym1      Asym2      xmid1      xmid2       scal 
 0.6035251  0.7129002 22.7821155 33.3080264  4.6928316

并且由于简化模型嵌套在完整模型中,我们可以再次进行似然比检验:

> anova(fm3, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym[group], xmid[group], scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1    115    0.13945                             
2    114    0.13934  1 0.00010637   0.087 0.7685

较大的p值表明缩减模型与完整模型一样符合预期。

我们当然可以做类似的测试来检查 just Asym、 justxmid或 both 是否需要不同的参数值。也就是说,我建议像这样进行逐步回归来消除参数。相反,只需fm2针对简单模型 ( ) 测试完整模型 ( fm1),并对结果感到满意。为了量化任何差异,绘图将很有帮助。