fit: Fit General Linear Model

In bhagwataditya/autonomics: Unified statistal Modeling of Omics Data

Fit General Linear Model

Description

Fit General Linear Model

Usage

fit(
  object,
  groupvar = "subgroup",
  formula = group2formula(groupvar, object),
  engine = "limma",
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = contr.treatment.explicit,
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun,
    verbose = FALSE),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) contrast_coefs(design = design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p", "se", "t")[1:2],
  ftest = if (is.null(coefs)) TRUE else FALSE,
  sep = FITSEP,
  suffix = paste0(sep, engine),
  verbose = TRUE,
  plot = FALSE
)

fit_limma(
  object,
  groupvar = "subgroup",
  formula = as.formula(sprintf("~ %s", groupvar)),
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = contr.treatment.explicit,
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) model_coefs(design = design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p", "t"),
  ftest = if (is.null(coefs)) TRUE else FALSE,
  sep = FITSEP,
  suffix = paste0(sep, "limma"),
  verbose = TRUE,
  plot = FALSE
)

.fit_limma(
  object,
  groupvar = "subgroup",
  formula = as.formula(sprintf("~ %s", groupvar)),
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = contr.treatment.explicit,
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) model_coefs(design = design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p", "se", "t")[1:2],
  ftest = if (is.null(coefs)) TRUE else FALSE,
  sep = FITSEP,
  suffix = paste0(sep, "limma"),
  verbose = TRUE,
  plot = FALSE
)

fit_wilcoxon(
  object,
  groupvar = "subgroup",
  formula = as.formula(sprintf("~ %s", groupvar)),
  drop = NULL,
  codingfun = contr.treatment.explicit,
  design = NULL,
  contrasts = NULL,
  coefs = NULL,
  block = NULL,
  weightvar = NULL,
  statvars = c("effect", "p"),
  sep = FITSEP,
  suffix = paste0(sep, "wilcoxon"),
  verbose = TRUE,
  plot = FALSE
)

Arguments

object	SummarizedExperiment
formula	model formula
engine	'limma', 'lm', 'lme', 'lmer', or 'wilcoxon'
drop	TRUE or FALSE
codingfun	factor coding function contr.treatment: intercept = y0, coefi = yi - y0 contr.treatment.explicit: intercept = y0, coefi = yi - y0 code_control: intercept = ymean, coefi = yi - y0 contr.diff: intercept = y0, coefi = yi - y(i-1) code_diff: intercept = ymean, coefi = yi - y(i-1) code_diff_forward: intercept = ymean, coefi = yi - y(i+) code_deviation: intercept = ymean, coefi = yi - ymean (drop last) code_deviation_first: intercept = ymean, coefi = yi - ymean (drop first) code_helmert: intercept = ymean, coefi = yi - mean(y0:(yi-1)) code_helmert_forward: intercept = ymean, coefi = yi - mean(y(i+1):yp)
design	design matrix
contrasts	NULL or character vector: coefficient contrasts to test
coefs	NULL or character vector: model coefs to test
block	block svar (or NULL)
weightvar	NULL or name of weight matrix in assays(object)
statvars	character vector: subset of c('effect', 'p', 'fdr', 't', 'se')
ftest	TRUE or FALSE
sep	string: pvar separator ("~" in "p~t2~limma")
suffix	string: pvar suffix ("limma" in "p~t2~limma")
verbose	whether to msg
plot	whether to plot

Value

Updated SummarizedExperiment

Examples

# Read
  file <- system.file('extdata/atkin.metabolon.xlsx', package = 'autonomics')
  object <- read_metabolon(file)
    
# Standard
  fdt(object) %<>% extract(, 'feature_id')
  object %<>% fit_lm(        ~ subgroup)                     #     statistics default
  object %<>% fit_limma(     ~ subgroup)                     # bioinformatics default
  summarize_fit(object)
    
# Blocked
  fdt(object) %<>% extract(, 'feature_id')
  object %<>% fit_limma(     ~ subgroup, block = 'Subject')  #        simple random effects
  object %<>% fit_lme(       ~ subgroup, block = 'Subject')  #      powerful random effects
  object %<>% fit_lmer(      ~ subgroup, block = 'Subject')  # more powerful random effects
  summarize_fit(object)
    
# Alternative coding: e.g. grand mean intercept
  fdt(object) %<>% extract(, 'feature_id')
  object %<>% fit_limma(     ~ subgroup, block = 'Subject', codingfun = code_control)
  object %<>% fit_lme(       ~ subgroup, block = 'Subject', codingfun = code_control)
  object %<>% fit_lmer(      ~ subgroup, block = 'Subject', codingfun = code_control)
  summarize_fit(object)
    
# Posthoc contrasts (only limma!)
  fdt(object) %<>% extract(, 'feature_id')
  object %<>% fit_limma( ~ subgroup, block = 'Subject', codingfun = code_control, coefs ='t1-t0')
  object %<>% fit_limma( ~ 0 + subgroup, block = 'Subject', contrasts = 't1-t0')
      # flexible, but only approximate
      # stat.ethz.ch/pipermail/bioconductor/2014-February/057682.html
        
# Non-parametric: wilcoxon
  fdt(object) %<>% extract(, 'feature_id')
  object %<>% fit_wilcoxon( ~ subgroup)                    # unpaired
  object %<>% fit_wilcoxon( ~ subgroup, block = 'Subject') # paired
    
# Custom separator
  fdt(object) %<>% extract(, 'feature_id')
  fdt( fit_lm(      object, sep = '.'))
  fdt( fit_limma(   object, block = 'Subject', sep = '.') )
  fdt( fit_lme(     object, block = 'Subject', sep = '.') )
  fdt( fit_lmer(    object, block = 'Subject', sep = '.') )
  fdt( fit_wilcoxon(object, block = 'Subject', sep = '.') )
  fdt( fit_wilcoxon(object, sep = '.') )

bhagwataditya/autonomics documentation built on Nov. 1, 2024, 5:47 a.m.