maxEM: maximum likelihood estimates relative abundances

In sonicLength: Estimating Abundance of Clones from DNA Fragmentation Data

Description

From information about the lengths of sonicants of integration sites, infer the relative abundances of different clones and the distributon of sonicant lengths

Usage

maxEM(slmat, theta.var=FALSE, phi.update=NULL,
                  phi.deriv=NULL, lframe = NULL, glm.frm = NULL, iter.control=NULL, ... ) 

Arguments

slmat	a matrix whose rows correspond to unique lengths with rownames indicating those lengths
theta.var	logical, return variance of theta estimates?
phi.update	a function of an object like slmat that returns estimtes of phi - a default vesrion is invoked if none is provided
phi.deriv	function of theta and phi that returns derivatives of phi wrt beta (its parameters)
lframe	a data.frame which will be used to estimate phi, the supplied function phi.update must know what to do with this arg or it will be ignored
glm.frm	a formula like y ~ bs(x,knots=c(50,100)) to fit phi
iter.control	a list of default values for iteration control - see maxEM.iter.control
...	possibly other args to pass to phi.update

Details

The EM method is used to infer the relative abundances of different sites.

Value

theta	a vector of the abundance estimates
phi	a vector of the probabilities of sonicant lengths
call	the call used

Author(s)

Charles C. Berry [e-mail]

See Also

Ey.given.x

Examples

mat0 <- matrix(0,nr=48,nc=140)
vals <- c(rep(1,100),2:40,100)
mat1 <- sapply( vals, function(x) as.numeric(is.element(1:200 ,rgeom(x,.02))))
mat <- rbind(mat0,mat1)
posVals <- colSums(mat) > 0
vals <- vals[ posVals ]
mat <- mat[, posVals ]
rownames(mat) <- 1:nrow(mat)
res <- maxEM(mat)
matplot( vals, cbind(res$theta, colSums(mat)), pch=1:2,
        xlab='original values', ylab='estimated values',
        main='Simulated Sonicants and Estimates')
legend( "bottomright", pch=1:2, col=1:2,
        legend=c(expression(hat(theta)[j]),expression(sum(y[ij],i))))
abline(a=0,b=1,col='gray')

sonicLength documentation built on Sept. 20, 2021, 9:06 a.m.