Hi all,
I'm trying to analyze a dataset that has left-censored, right-censored, and
categorical predictor variables and a univariate response in mcmcglmm. However,
I am uncertain of how to specify the appropriate model command and priors.
The response (rates) is a net diversification rate derived from BAMM. The
right-censored data is body size (max.min, max=Inf), and size at sexual
maturity is left censored (f.min= -Inf, f.max). There are four binary variables
indicating presence/absence: M.sup, M.inf, F.sup, F.inf.
My code is as follows:
library(MCMCglmm)
dat <- read.table("appended3.txt", row.names=1, header=TRUE) ##Read data table
dat$f.min[which(is.na(dat$f.min))]<- -Inf ##Assign -Inf for left-censored
information
t <- read.tree("MCC_pruned.tre") ##Read phylogenetic tree
tnames <- t$tip.label[order(t$tip.label)] ##Get tip names
full <- as.data.frame(cbind(tnames, dat)) ##Bind tip names to existing data
frame
invT <- inverseA(t)$Ainv ##Calculate inverse, assuming nodes="ALL", the
default
prior <- list(R=list(V=1,nu=0.002),
G=list(G1=list(V=1,
nu=0.002))) ##General
"get it running" prior I have seen in tutorials
model_test <- MCMCglmm(rates ~ cbind(max.min, max) + cbind(f.min, f.max)+ M.sup
+ M.inf + F.sup + F.inf,
random = ~tnames, prior=prior,
family=c("gaussian","cenexponential","cenexponential", "threshold","threshold",
"threshold","threshold"), ginverse=list(tnames=invT),
data=full, nitt=5000000,burnin=1000,thin=500)
##The above specification is attempting to assign gaussian for rates,
cenexponential for the left-and right-censored data, and threshold for the four
binary characters.
I get the error
Error in matrix(unlist(value, recursive = FALSE, use.names = FALSE), nrow = nr,
:
'data' must be of a vector type, was 'NULL'
traceback() gives the following:
3: matrix(unlist(value, recursive = FALSE, use.names = FALSE), nrow = nr,
dimnames = list(rn, cn))
2: Ops.data.frame(data[, which(names(data) == response.names[nt +
1])], data[, which(names(data) == response.names[nt])])
1: MCMCglmm(rates ~ cbind(max.min, max) + cbind(f.min, f.max) +
M.sup + M.inf + F.sup + F.inf, random = ~tnames, prior = prior,
family = c("cenexponential", "cenexponential", "threshold",
"threshold", "threshold", "threshold"), ginverse = list(tnames =
invT),
data = full, nitt = 5e+06, burnin = 1000, thin = 500)
and debug() stops after
debug: if (any(data[, which(names(data) == response.names[nt + 1])] <
data[, which(names(data) == response.names[nt])], na.rm = T)) {
stop("for censored traits left censoring point must be less than right
censoring point")
}
However, the lower bound for max is equal to or larger than the upper bound for
f.min.
How can I get this running?
Many thanks,
~John
John S. S. Denton, Ph.D.
Department of Vertebrate Paleontology
American Museum of Natural History
www.johnssdenton.com
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
