0

After fitting a non-linear mixed model with crossed categorical covariates, how do you write the model equation (with parameter estimates) for each group combination?

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

data(CO2)

Using the CO2 data set as referenced in Pinheiro and Bates 2000:

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

Four-paneled figure, breaking down by Type and Treatment, plotting uptake against conc

Fit a model with an Asymptotic regression model with an offset, allowing the Asym and lrc to vary by Type and Treatment (but requiring c0 be constant), according to a model in section 8.2.2 in Pinheiro and Bates.

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

(output truncated)

I know how to plot the modeled data with conditional predicted values (using the plant information "Plant") and marginal predicted values (ignoring the plant information; "fixed")

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

enter image description here

And I know how to add onto data frame the marginal predicted values (ignoring the effect of the plant random effect)

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

Is there an easy way to obtain the parameters for the four model equations for each of the Type * Treatment combos?

emudrak
  • 789
  • 8
  • 25

1 Answers1

0

Yes! Thanks the to always-useful emmeans package, this is quite straightforward! For each parameter, make sure the model equation after the ~ matches that of the fixed list of the original model fit.

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

 Type        Treatment  emmean   SE df lower.CL upper.CL
 Quebec      nonchilled   41.8 1.48 64     38.9     44.8
 Mississippi nonchilled   31.3 1.48 64     28.3     34.2
 Quebec      chilled      38.8 1.49 64     35.9     41.8
 Mississippi chilled      17.4 1.44 64     14.5     20.3


Confidence level used: 0.95

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

 Type        Treatment  emmean     SE df lower.CL upper.CL
 Quebec      nonchilled  -4.56 0.0910 64    -4.74    -4.38
 Mississippi nonchilled  -4.66 0.1015 64    -4.86    -4.46
 Quebec      chilled     -4.73 0.0897 64    -4.91    -4.55
 Mississippi chilled     -4.09 0.1715 64    -4.43    -3.75

Confidence level used: 0.95

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

 1       emmean   SE df lower.CL upper.CL
 overall   50.5 4.12 64     42.3     58.7

Confidence level used: 0.95

Gather these values in a single data frame:

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
emudrak
  • 789
  • 8
  • 25