maxEM: maximum likelihood estimates relative abundances

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

Value

theta

a vector of the abundance estimates

phi

a vector of the probabilities of sonicant lengths

call

the call used

Details

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

Examples

Run this code

# NOT RUN {
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')
# }