具有范畴协变量的非线性混合模型:按群求参数估计

Question

在拟合具有交叉范畴协变量的非线性混合模型后，如何为每一组组合编写模型方程(带有参数估计)？

library(ggplot2)
library(tidyr)
library(dplyr)

data(CO2)

使用CO2数据集作为皮涅罗和贝茨2000年中的引用

ggplot(CO2, aes(x=conc, y=uptake, group=Plant))+
  geom_line()+
  facet_grid(Treatment ~ Type)

​

​

根据皮涅罗和贝茨8.2.2节中的模型，用带有偏移量的渐近回归模型拟合模型，允许Asym和lrc随Type和Treatment的变化而变化(但要求c0是常数)。

library(nlme)
fm4CO2.nlme <- nlme(uptake ~ SSasympOff(conc, Asym, lrc, c0), data=CO2, 
                    fixed= list(Asym + lrc ~ Type * Treatment, c0 ~ 1),
                    random = Asym + lrc ~ 1,
                    start = c(32.412, 0, 0, 0, -4.5603,0, 0, 0, 49.344))

summary(fm4CO2.nlme)

Nonlinear mixed-effects model fit by maximum likelihood
  Model: uptake ~ SSasympOff(conc, Asym, lrc, c0) 
  Data: CO2 
      AIC      BIC   logLik
  388.418 420.0186 -181.209

Random effects:
 Formula: list(Asym ~ 1, lrc ~ 1)
 Level: Plant
 Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr  
Asym.(Intercept) 2.34968017 As.(I)
lrc.(Intercept)  0.07960176 -0.92 
Residual         1.79195963       

Fixed effects:  list(Asym + lrc ~ Type * Treatment, c0 ~ 1) 
                                          Value Std.Error DF   t-value p-value
Asym.(Intercept)                       41.81755  1.562449 64  26.76410  0.0000
Asym.TypeMississippi                  -10.53047  2.208351 64  -4.76848  0.0000
Asym.Treatmentchilled                  -2.96942  2.213205 64  -1.34169  0.1844
Asym.TypeMississippi:Treatmentchilled -10.89926  3.112279 64  -3.50202  0.0008
lrc.(Intercept)                        -4.55726  0.096292 64 -47.32762  0.0000
lrc.TypeMississippi                    -0.10411  0.121685 64  -0.85557  0.3954
lrc.Treatmentchilled                   -0.17124  0.111962 64  -1.52947  0.1311
lrc.TypeMississippi:Treatmentchilled    0.74124  0.221716 64   3.34322  0.0014
c0                                     50.50804  4.364848 64  11.57155  0.000

(输出截断)

我知道如何用条件预测值(使用植物信息“植物”)和边际预测值(忽略植物信息；“固定”)绘制建模数据。

plot(augPred(fm4CO2.nlme, level=0:1))

​

​

我知道如何在数据框架中添加边际预测值(忽略植物随机效应的影响)。

CO2$predict=predict(fm4CO2.nlme, CO2 , level=0)

对于每一个Type * Treatment 组合，是否有一种简单的方法来获得四个模型方程的参数？

Answer 1

是!感谢这个总是有用的emmeans包，这是非常简单的！对于每个参数，在~与原始模型的fixed列表匹配之后，确保模型方程匹配。

library(emmeans)
emmeans(fm4CO2.nlme,  ~ Type * Treatment, param="Asym")

lower.CL upper.CL魁北克非冰镇41.8 1.48 64 38.9 44.8密西西比州非冰鲜31.3 1.48 28.3 34.2魁北克冰镇38.8 1.49 64 35.9 41.8密西西比州冷冻17.4 1.44 64 14.5 20.3置信水平: 0.95

 emmeans(fm4CO2.nlme,  ~ Type * Treatment, param="lrc")

lower.CL upper.CL魁北克省非冰镇-4.56 0.0910 64 -4.74 -4.38密西西比非冰镇-4.66 0.1015 64 -4.86 -4.46魁北克-4.73 64 -4.91 -4.55密西西比州冰镇-4.09 0.1715 64 -4.43 -3.75置信度: 0.95

emmeans(fm4CO2.nlme,  ~ 1, param="c0")

1平均SE df lower.CL upper.CL总体使用50.5 4.12 64 42.3 58.7置信度: 0.95

在一个数据框架中收集这些值：

Asym_vals <- emmeans(fm4CO2.nlme,  ~ Type * Treatment, param="Asym") %>%
  as.data.frame()%>%
  rename_with(.fn=~paste0(., ".Asym"), .cols=emmean:upper.CL)

lrc_vals <- emmeans(fm4CO2.nlme,  ~ Type * Treatment, param="lrc") %>%
  as.data.frame()%>%
  rename_with(.fn=~paste0(., ".lrc"), .cols=emmean:upper.CL)

c0_vals <- emmeans(fm4CO2.nlme,  ~ 1, param="c0") %>%
  as.data.frame()%>%
  rename_with(.fn=~paste0(., ".c0"), .cols=emmean:upper.CL)


modeleqns <- Asym_vals %>% 
  left_join(lrc_vals, by=c("Type", "Treatment"), suffix=c(".Asym", ".lrc"))%>%
  merge(c0_vals, by=NULL)%>%
  select(Type, Treatment, starts_with("emmean."))
modeleqns

Type Treatment emmean.Asym emmean.lrc emmean.c0 1 Quebec nonchilled 41.81755 -4.557257 50.50804 2 Mississippi nonchilled 31.28708 -4.661367 50.50804 3 Quebec chilled 38.84812 -4.728499 50.50804 4 Mississippi chilled 17.41839 -4.091365 50.50804