negloglik_adjusted_profile: negative adjusted profile log-likelihood function of a CARseq...
In Sun-lab/CARseq: Cell type-specific differential analysis of count data

Description Usage Arguments References Examples

negloglik_adjusted_profile provides a negative adjusted profile log-likelihood function of negative binomial distribution whose mean is a sum of nonnegative terms with covariates. This is used to estimate the overdispersion parameter.

negloglik_adjusted_profile(
  logtheta,
  mu_matrix,
  coefs,
  cell_type_specific_variables,
  other_variables,
  read_depth,
  cellular_proportions,
  counts,
  is_active = NULL,
  use_log_scale_params = FALSE
)

`logtheta`	log(overdispersion) as a scalar. Overdispersion is parametrized as `theta` in glm.nb.
`mu_matrix`	cell type-specific fitted values from CARseq negative binomial model.
`coefs`	a vector of c(beta, gamma, overdispersion), where beta is a length K vector of cell type-indepedent coefficients, gamma is a matrix of H x M dimension of cell type-specific non-negative coefficients, and overdispersion is a scalar. Note: only gamma is used in the evaluation here.
`cell_type_specific_variables`	an array of n_B x H x K of cell type-independent variables.
`other_variables`	a design matrix of n_B x M of cell type-specific variables. Can be NULL.
`read_depth`	a vector of n_B sample-specific read depths.
`cellular_proportions`	a matrix of n_B x H of cellular proportions.
`counts`	A vector of n_B total read counts observed.
`is_active`	A logical vector recording which predictors are included in the model.
`use_log_scale_params`	logical. The default is FALSE, where the cell type-specific effects are parametrized on a non-log scale. Otherwise, they are parametrized on a log scale.

Davis J. McCarthy, Yunshun Chen, Gordon K. Smyth, Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation, Nucleic Acids Research, Volume 40, Issue 10, 1 May 2012, Pages 4288-4297, https://doi.org/10.1093/nar/gks042

library(CARseq)
set.seed(1234)
H = 4
n_B = 60
K = 1
M = 3
overdispersion_theta = 10
coefs = c(rep(0, K), rep(1, H*M), overdispersion_theta)  # initial value
cell_type_specific_variables = array(0, dim=c(n_B, H, M))
cell_type_specific_variables[1:(n_B/3), , 1] = 1
cell_type_specific_variables[(n_B/3+1):(n_B*2/3), , 2] = 1
cell_type_specific_variables[(n_B*2/3+1):n_B, , 3] = 1
other_variables = matrix(runif(n_B * K, min = -1, max = 1), nrow=n_B, ncol=K)
read_depth = round(runif(n_B, min = 50, max = 100))
cellular_proportions = matrix(runif(n_B * H), nrow = n_B, ncol = H)
cellular_proportions = cellular_proportions / rowSums(cellular_proportions)
counts = round(runif(n_B, min = 100, max = 200))
is_active = rep(TRUE, K + H*M)
res = CARseq:::negloglik(coefs, cell_type_specific_variables, other_variables, read_depth, cellular_proportions, counts, verbose=TRUE)
CARseq:::negloglik_adjusted_profile(log(overdispersion_theta), attr(res, "cell.type.specific.fitted.values"), coefs, cell_type_specific_variables, other_variables, read_depth, cellular_proportions, counts, is_active)