Nothing
## this is the heart of model selection and
## profile likelihood
pva.llr <-
function(null, alt, pred)
{
if (!inherits(null, "pva"))
stop("'null' must be of class 'pva'")
if (!inherits(alt, "pva"))
stop("'alt' must be of class 'pva'")
err0 <- null@model@obs.error
err1 <- alt@model@obs.error
obs <- null@observations
if (!identical(obs, alt@observations))
stop("data in null amd alternative model must be identical")
if (err0 != "none" && err1 == "none")
stop("switch null and alternative model")
if (identical(null@model@growth.model, alt@model@growth.model) &&
identical(null@model@obs.error, alt@model@obs.error) &&
(length(null@coef) > length(alt@coef)))
warning("Hint: check what is null and alternative ",
"if models are nested")
if (!any(is.na(obs)) && err1 == "none" && missing(pred)) {
pred <- matrix(log(obs), nrow=1)
}
data0 <- switch(err0,
"none" = log(obs),
"normal" = log(obs),
"poisson" = obs)
data1 <- switch(err1,
"none" = log(obs),
"normal" = log(obs),
"poisson" = obs)
## use here parApply with parallel package -- future stuff
logd0 <- apply(pred, 1, null@model@logdensity,
mle=coef(null), data=data0,
null_obserror=err0 != "none", alt_obserror=err1 != "none")
logd1 <- apply(pred, 1, alt@model@logdensity,
mle=coef(alt), data=data1,
null_obserror=err0 != "none", alt_obserror=err1 != "none")
## log likelihood ratio
if (err0 == "none" && err1 != "none" && any(is.na(obs))) {
# out <- log((1/(mean(1/exp(logd0)))) / (1/(mean(1/exp(logd1))))) ## --------------check this
stop("missing values not yet fully supported in model comparison")
} else {
out <- log(mean(exp(logd0 - logd1)))
}
out
}
## model selection
model.select <-
function(null, alt, B=10^4)
{
## generateLatent takes care of the non-obs.err and
## missing value cases
pred <- suppressWarnings(generateLatent(alt,
n.chains=1, n.iter=B))
llr <- pva.llr(null, alt, pred)
## neffective can be different (e.g. Ricer vs. Gompertz)
n0 <- null@model@neffective(null@observations)
n1 <- alt@model@neffective(null@observations)
p0 <- null@model@p - length(null@model@fixed)
p1 <- alt@model@p - length(alt@model@fixed)
rval <- data.frame(log_LR = llr,
delta_AIC = -2*llr + 2*(p0-p1),
delta_BIC = -2*llr + log(n0)*p0 - log(n1)*p1,
delta_AICc= -2*llr + 2*((p0-p1) +
p0*(p0+1)/(n0-p0-1) - p1*(p1+1)/(n0-p1-1)))
attr(rval, "fancy") <- list(
data=fancyPVAmodel(null, "", part=2),
null=fancyPVAmodel(null, "", part=1),
alt=fancyPVAmodel(alt, "", part=1))
attr(rval, "models") <- list(
null = list(name=substitute(null), n=n0, p=p0,
growth.model=null@model@growth.model,
obs.error = null@model@obs.error),
alt = list(name=substitute(alt), n=n1, p=p1,
growth.model=alt@model@growth.model,
obs.error = alt@model@obs.error))
class(rval) <- c("pvaModelSelect", "data.frame")
rval
}
print.pvaModelSelect <-
function(x, ...)
{
f <- attr(x, "fancy")
m <- attr(x, "models")
cat("PVA Model Selection:\n")
cat(f$data, "\n\n")
cat("Null Model:", m$null$name, "\n ", f$null, "\n\n")
cat("Alternative Model:", m$alt$name, "\n ", f$alt, "\n\n")
z <- x
class(z) <- "data.frame"
print(z, ...)
llr <- x$log_LR
if (abs(llr) < 10^-6) {
cat("\nNull and Alternative models are",
"equally well supported\n\n")
} else {
if (-llr < 0) {
good <- "Null Model"
bad <- "Alternative Model"
} else {
good <- "Alternative Model"
bad <- "Null Model"
}
if (abs(llr) < 2) {
qualifier <- "slightly better"
} else if (abs(llr) < 8) {
qualifier <- "better"
} else qualifier <- "strongly"
cat("\n", good, " is ", qualifier, " supported over the ",
bad, "\n\n", sep="")
}
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.