run_reduced_em_algo: Run univariate poisson EM algo

In timothy-barry/glmeiv: Implements the GLM-EIV model and method

Run univariate poisson EM algo

Description

This function aims to estimate three parameters: m_perturbation, g_perturbation, and pi.

Usage

run_reduced_em_algo(
  m,
  g,
  m_fitted,
  g_fitted,
  m_pert_guess,
  g_pert_guess,
  pi_guess,
  m_fam,
  g_fam,
  ep_tol = 0.5 * 1e-04,
  max_it = 50,
  min_it = 3,
  use_mrna_modality = TRUE
)

Arguments

m	mRNA counts
g	gRNA counts
m_fitted	fitted values on linear scale of mRNA model
g_fitted	fitted values on linear scale of gRNA model
m_pert_guess	initial guess for m perturbation parameter
g_pert_guess	initial guess for g perturbation parameter
pi_guess	initial guess for pi
m_fam	family object describing mRNA distribution
g_fam	family object describing gRNA distribution
ep_tol	relative tolerance limit for EM algorithm
max_it	maximum number of iterations before giving up
min_it	mininum number of iterations before declaring convergence

Details

The function takes offsets for the mRNA and gRNA modalities. It is assumed that all auxilliary information (baseline expression level, library size, batch effect and other technical factors, etc.) has been distilled into the offset vectors.

The function requires initial estimates for the parameters m_perturbation, g_perturbation, and pi.

Value

a list containing (i) the vector of estimates, (ii) log-likelihood of the fitted model, and (iii) the posterior membership probabilities.

Examples

set.seed(4)
# NB response
m_fam <- g_fam <- augment_family_object(MASS::negative.binomial(20))
m_intercept <- 0
g_intercept <- 0
m_perturbation <- log(0.6)
g_perturbation <- log(1.4)
pi <- 0.02
n <- 100000
m_offset <- log(stats::rpois(n, 100))
g_offset <- log(stats::rpois(n, 50))
dat <- generate_full_data(m_fam = m_fam, m_intercept = m_intercept, m_perturbation = m_perturbation,
g_fam = g_fam, g_intercept = g_intercept, g_perturbation = g_perturbation, pi = pi, n = n, B = 2,
covariate_matrix = NULL, m_covariate_coefs = NULL, g_covariate_coefs = NULL, m_offset = m_offset,
g_offset = g_offset)[[1]]
m_fit <- glm(formula = m ~ p + 0, family = m_fam, data = dat, offset = m_offset)
g_fit <- glm(formula = g ~ p + 0, family = g_fam, data = dat, offset = g_offset)
m_pert_guess <- log(runif(1, 0.1, 1.5))
g_pert_guess <- log(runif(1, 0.5, 10))
pi_guess <- runif(1, 0, 0.05)
m <- dat$m
g <- dat$g
m_fitted <- m_offset
g_fitted <- g_offset
fit_univariate <- run_reduced_em_algo(m, g, m_fitted, g_fitted,
m_pert_guess, g_pert_guess, pi_guess, m_fam, g_fam)

# Gaussian
m_fam <- g_fam <- augment_family_object(gaussian())
m_intercept <- 0
g_intercept <- 0
m_perturbation <- 2
g_perturbation <- -1
pi <- 0.02
n <- 100000
m_offset <- log(stats::rpois(n, 100))
g_offset <- log(stats::rpois(n, 50))
dat <- generate_full_data(m_fam = m_fam, m_intercept = m_intercept, m_perturbation = m_perturbation,
g_fam = g_fam, g_intercept = g_intercept, g_perturbation = g_perturbation, pi = pi, n = n, B = 2,
covariate_matrix = NULL, m_covariate_coefs = NULL, g_covariate_coefs = NULL, m_offset = m_offset,
g_offset = g_offset)[[1]]
m_fit <- glm(formula = m ~ p + 0, family = m_fam, data = dat, offset = m_offset)
g_fit <- glm(formula = g ~ p + 0, family = g_fam, data = dat, offset = g_offset)
m_pert_guess <- 3
g_pert_guess <- -2
pi_guess <- runif(1, 0, 0.05)
m <- dat$m
g <- dat$g
m_fitted <- m_offset
g_fitted <- g_offset
fit_univariate <- run_reduced_em_algo(m, g, m_fitted, g_fitted,
m_pert_guess, g_pert_guess, pi_guess, m_fam, g_fam)

timothy-barry/glmeiv documentation built on Jan. 30, 2024, 3:46 p.m.