EstimateCCM: Estimate the parameters in CCM

In beiko-lab/evolCCM: R package for Community Coevolution Model

Estimate the parameters in CCM

Description

The function for estimating the parameters in CCM.

Usage

EstimateCCM = function(profiles, phytree, ip=0.1, pen=0.5, ... )

Arguments

profiles	a matrix containing the profiles. Columns are profiles and rows are species.
phytree	a phylogenetic tree.
ip	the initial values for optimizer. For a better convergence, a good set of starting values could be the estimates by a large tuning parameter.
pen	the tuning parameter \lambda. Default value is 0.5. Value of 0 means no regularization.
...	control parameters available to optimizer ‘nlminb' such as 'trace', 'rel.tol', .... .Example See '?nlminb', ’control' argument .

Value

a list with following elements

• alpha: estimated intrinsic rates.

• B: estimated association matrix.

• nlm.par: all estimated parameters in the order as 'c(alpha, diag(B), B[upper.tri(B)])'.

• nlm.converge: convergence message. See '?nlminb'.

• nlm.hessian: estimated Hessian matrix used for calculating standard errors.

References

Paradis E, Schliep K (2019). “ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R.” Bioinformatics, 35, 526-528.

David M. Gay (1990), Usage summary for selected optimization routines. Computing Science Technical Report 153, AT&T Bell Laboratories, Murray Hill.

Examples

set.seed(123)
# generate a ranom 200-tip tree.
# Larger tree contains more information and tends to give better MLEs.
t <- rtree(200)
# setting arbitrary underlying parameters for a pair
n <- 2
alpha <- c(0.1, 0.1)
B <- matrix(0, n, n)
B[1,2] <- B[2,1] <- 1
diag(B) <- c(-0.3, 0.3)

# using the same set of parameter to simulate 20 pairs
# and estimate with CCM.

trueP = c(alpha, diag(B), B[upper.tri(B)]) # true parameters
estP = matrix(NA, nrow=20, ncol=length(trueP))
for (i in 1:20){
    simDF <- SimulateProfiles(t, alpha, B)
    aE <- EstimateCCM(profiles = simDF, phytree=t)
    estP[i,] = c(aE$alpha, diag(aE$B), aE$B[upper.tri(aE$B)])
    # or estP[i,] = aE$par same order as above
    print(i)
}
# plot the distribution of estimates
boxplot(estP)
points(1:length(trueP), trueP, pch=8, col="red")

beiko-lab/evolCCM documentation built on Feb. 26, 2024, 5:21 p.m.