npem.em: Fit the normal-Poisson model to data on a cell proliferation...

In kbroman/npem: Analysis of Cell Proliferation Assays

Fit the normal-Poisson model to data on a cell proliferation assay

Description

Uses a version of the EM algorithm to fit the normal-Poisson mixture model to data on a cell proliferation assay.

Usage

npem.em(
  y,
  ests,
  cells = 10^6,
  n = c(24, 24, 24, 22),
  n.plates = 1,
  use.order.constraint = TRUE,
  maxit = 2000,
  tol = 0.000001,
  maxk = 20,
  prnt = 0
)

Arguments

y	Vector of transformed scintillation counts, in lexicographical order (plate by plate and group by group within a plate.)
ests	Initial parameter estimates, as a vector of length n.groups + 3*n.plates, of the form (\lambda's, (a, b, \sigma)'s), where \lambda is the average number of responding cells per 10^6 cells for a group, and (a, b, \sigma) are the plate-specific parameters.
cells	Number of cells per well. The \lambda's will be rescaled to give response per 10^6 cells. This may be either a single number (if all wells have the same number of cells, or 10^6 if one wishes the \lambda's to not be rescaled), a value for each plate (vector of length n.plates, or a value for each well (a vector of the same length as y).
n	Vector giving the number of wells within each group. This may have length either n.groups (if all plates have the same number of wells per group) or n.groups*n.plates.
n.plates	The number of plates in the data.
use.order.constraint	If TRUE, force the constraint \lambda_0 \le \lambda_i for all i \ge 1; otherwise, no constraints are applied.
maxit	Maximum number of EM iterations to perform.
tol	Tolerance to determine when to stop the EM algorithm.
maxk	Maximum k value in sum calculating E(k | y).
prnt	If 0, don't print anything; if 1, print the log likelihood at each iteration; and if 2, print the est's and the log lik. at each iteration.

Details

Calculations are performed in a C routine. [I should describe the normal-Poisson mixture model here.]

Value

ests	The estimated parameters in same form as the input argument ests.
k	The estimated number of responding cells per well, E(k|y).
ksq	E(k^2|y)
loglik	The value of the log likelihood at ests.
n.iter	Number of iterations performed.

Author(s)

Karl W Broman, [e-mail]

References

Broman et al. (1996) Estimation of antigen-responsive T cell frequencies in PBMC from human subjects. J Immunol Meth 198:119-132
Dempster et al. (1977) Maximum likelihood estimation from incomplete data via the EM algorithm. J Roy Statist Soc Ser B 39:1-38

See Also

npem.sem(), npem.start(), npsim(), npem.ll()

Examples

  # get access to an example data set
  data(p713)

  # analysis of plate3
    # get starting values
  start.pl3 <- npem.start(p713$counts[[3]],n=p713$n)
    # get estimates
  out.pl3 <- npem.em(p713$counts[[3]],start.pl3,n=p713$n)
    # look at log likelihood at starting and ending points
  npem.ll(p713$counts[[3]],start.pl3,n=p713$n)
  npem.ll(p713$counts[[3]],out.pl3$ests,n=p713$n)
    # repeat with great precision, starting at previous endpoint
  out.pl3 <- npem.em(p713$counts[[3]],out.pl3$ests,
                     n=p713$n,tol=1e-13)
    # run SEM algorithm to get standard errors
  out.sem.pl3 <- npem.sem(p713$counts[[3]],out.pl3,n=p713$n)
  round(out.pl3$ests,3)
  round(out.sem.pl3$se,3)

  # repeat the above for the pair, plates 3 and 4
    # get starting values
  start.pl34 <- npem.start(unlist(p713$counts[3:4]),n=p713$n,n.plates=2)
    # get estimates
  out.pl34 <- npem.em(unlist(p713$counts[3:4]),start.pl34,n=p713$n,n.plates=2)
    # look at log likelihood at starting and ending points
  npem.ll(unlist(p713$counts[3:4]),start.pl34,n=p713$n,n.plates=2)
  npem.ll(unlist(p713$counts[3:4]),out.pl34$ests,n=p713$n,n.plates=2)
    # repeat with great precision, starting at previous endpoint
  out.pl34 <- npem.em(unlist(p713$counts[3:4]),out.pl34$ests,
                     n=p713$n,tol=1e-13,n.plates=2)
    # run SEM algorithm to get standard errors
  out.sem.pl34 <- npem.sem(unlist(p713$counts[3:4]),out.pl34,n=p713$n,n.plates=2)
  round(out.pl34$ests,3)
  round(out.sem.pl34$se,3)

kbroman/npem documentation built on May 17, 2023, 11:52 p.m.