Refit mixed effects model

Based on a 'GlmmSeq' or 'lmmSeq' class result object, this function attempts to refit an identical model for a specific gene based on the data and fitting parameters stored in the results object and refitting using either lme4::glmer() for GlmmSeq objects or lmer() for lmmSeq objects. The fitted model can then be passed on to other packages such as emmeans to look at estimated marginal means for the model.

Usage

glmmRefit(object, gene, ...)

# S3 method for lmmSeq
glmmRefit(object, gene, formula = object@formula, ...)

# S3 method for GlmmSeq
glmmRefit(
  object,
  gene,
  formula = object@formula,
  control = object@info$control,
  family = NULL,
  ...
)

Arguments

object

A fitted results object of class GlmmSeq or lmmSeq

gene

A character value specifying a single gene to extract a fitted model for

...

Optional arguments passed to either lme4::glmer() or lme4::lmer()

formula

Optional formula to use when refitting model

control

Optional control parameters, see lme4::lmerControl() or lme4::glmerControl()

family

Optional GLM family when refitting GLMM using lme4::glmer() or glmmTMB()

Value

Fitted model of class lmerMod in the case of LMM, or glmerMod or glmmTMB for a GLMM dependent on the original method.

Examples

data(PEAC_minimal_load)
disp <- apply(tpm, 1, function(x) {
  (var(x, na.rm = TRUE)-mean(x, na.rm = TRUE))/(mean(x, na.rm = TRUE)**2)
})
glmmtest <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                     countdata = tpm[1:2, ],
                     metadata = metadata,
                     dispersion = disp,
                     verbose = FALSE)

# show summary for single gene
summary(glmmtest, "MS4A1")
#> Generalised linear mixed model
#> Method: lme4
#> Family: Negative Binomial
#> Formula: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
#> <environment: 0x7fd12830a4e0>
#>  Dispersion         AIC      logLik     meanExp 
#>    3.789428  898.682864 -441.341432    2.591740 
#> 
#> Fixed effects:
#>                                Estimate Std. Error
#> (Intercept)                     2.73637    0.33669
#> Timepoint                       0.04908    0.08941
#> EULAR_6mModerate                0.19330    0.51103
#> EULAR_6mNon-response           -1.72376    0.67780
#> Timepoint:EULAR_6mModerate     -0.21080    0.13478
#> Timepoint:EULAR_6mNon-response  0.17462    0.16558
#> 
#> Analysis of Deviance Table (Type II Wald chisquare tests)
#>                      Chisq Df Pr(>Chisq)
#> Timepoint          0.01125  1    0.91554
#> EULAR_6m           5.68288  2    0.05834
#> Timepoint:EULAR_6m 5.43787  2    0.06594

# refit a single model using lme4::glmer()
fit <- glmmRefit(glmmtest, "MS4A1")
#> boundary (singular) fit: see help('isSingular')

# refit model with reduced formula
fit2 <- glmmRefit(glmmtest, "MS4A1",
                  formula = count ~ Timepoint + EULAR_6m + (1 | PATID))

# LRT
anova(fit, fit2)
#> Data: data
#> Models:
#> fit2: count ~ Timepoint + EULAR_6m + (1 | PATID)
#> fit: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
#>      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#> fit2    6 899.53 916.41 -443.77   887.53                       
#> fit     8 898.68 921.18 -441.34   882.68 4.8505  2    0.08846 .
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1