setClass("unmarkedFitList",
representation(fits = "list"),
validity = function(object) {
fl <- object@fits
testY <- function(fit) {
f <- fit@formula
umf <- getData(fit)
D <- unmarked:::getDesign(umf, f)
D$y
}
umf1 <- getData(fl[[1]])
form1 <- fl[[1]]@formula
y1 <- unmarked:::getDesign(umf1, form1)$y
dataTest <- sapply(fl, function(x) isTRUE(all.equal(umf1, getData(x))))
yTest <- sapply(fl, function(x) isTRUE(all.equal(y1, testY(x))))
if(!all(dataTest)) {
stop("Data are not the same among models. Make sure you use the same unmarkedFrame object for all models.")
}
else if(!all(yTest)) {
stop("Data are not the same among models due to missing covariate values. Consider removing NAs before analysis.")
}
TRUE
}
)
# constructor of unmarkedFitList objects
fitList <- function(..., fits) {
if(length(list(...)) > 0 & !missing(fits))
stop("Do not use both the '...' and 'fits' arguments")
if(missing(fits)) {
fits <- list(...)
isList <- sapply(fits, function(x) is.list(x))
if(sum(isList) > 1)
stop("Specify models as common-seperated objects, or use fits = 'mylist'")
if(isList[1L]) {
warning("If supplying a list of fits, use fits = 'mylist'")
fits <- fits[[1L]] # This is allowed for back-compatability.
}
if(is.null(names(fits))) {
c <- match.call(expand.dots = FALSE)
names(fits) <- as.character(c[[2]])
warning("Your list was unnamed, so model names were added as object names")
}
}
if(is.null(names(fits))) {
names(fits) <- as.character(1:(length(fits)))
warning("Your list was unnamed, so model names were added as c('1','2',...)")
}
umfl <- new("unmarkedFitList", fits=fits)
return(umfl)
}
setMethod("summary", "unmarkedFitList", function(object) {
fits <- object@fits
for(i in 1:length(fits))
summary(fits[[i]])
})
setMethod("predict", "unmarkedFitList", function(object, type, newdata=NULL,
backTransform = TRUE, appendData = FALSE) {
fitList <- object@fits
ese <- lapply(fitList, predict, type = type, newdata = newdata,
backTransform = backTransform)
E <- sapply(ese, function(x) x[,"Predicted"])
SE <- sapply(ese, function(x) x[,"SE"])
ic <- sapply(fitList, slot, "AIC")
deltaic <- ic - min(ic)
wts <- exp(-deltaic / 2)
wts <- wts / sum(wts)
parav <- as.numeric(E %*% wts)
seav <- as.numeric((SE + (E - parav)^2) %*% wts) # Double check this
out <- data.frame(Predicted = parav, SE = seav)
if(appendData) {
if(missing(newdata))
newdata <- getData(object@fits[[1]])
out <- data.frame(out, newdata)
}
return(out)
})
# Condition number
cn <- function(object) {
h <- hessian(object)
if(any(is.na(h))) return(NA)
else {
ev <- eigen(h)$value
return(max(ev) / min(ev))
}
}
# R-squared index from Nagelkerke (1991)
nagR2 <- function(fit, nullfit)
{
n <- sampleSize(fit)
devI <- 2 * fit@negLogLike
devN <- 2 * nullfit@negLogLike
r2 <- 1 - exp((devI - devN) / n)
r2max <- 1 - exp(-1 * devN / n)
return(r2 / r2max)
}
setGeneric("modSel",
def = function(object, ...) {
standardGeneric("modSel")
}
)
setClass("unmarkedModSel",
representation(
Full = "data.frame",
Names = "matrix"
)
)
# Model selection results from an unmarkedFitList
setMethod("modSel", "unmarkedFitList",
function(object, nullmod=NULL)
{
if (!is.character(nullmod) && !is.null(nullmod)) {
stop("nullmod must be character name of null model fit in the fitlist.")
}
fits <- object@fits
estList <- lapply(fits, coef, altNames=TRUE)
seList <- lapply(fits, function(x)
if(any(is.na(x@opt$hessian))) rep(NA, length(coef(x)))
else sqrt(diag(vcov(x, altNames=TRUE))))
eNames <- sort(unique(unlist(sapply(estList, names))))
seNames <- paste("SE", eNames, sep="")
eseNames <- character(l <- length(c(eNames, seNames)))
eseNames[seq(1, l, by=2)] <- eNames
eseNames[seq(2, l, by=2)] <- seNames
cNames <- c("model", "formula", eseNames)
out <- data.frame(matrix(NA, ncol=length(cNames), nrow=length(fits)))
colnames(out) <- cNames
out$model <- names(fits)
out$formula <- sapply(fits, function(x) {
f <- as.character(x@formula)
f <- paste(f[2], "~", f[3])
f
})
for(i in 1:length(eNames)) {
out[,eNames[i]] <- sapply(estList, function(x) x[eNames[i]])
out[,seNames[i]] <- sapply(seList, function(x) x[eNames[i]])
}
out$Converge <- sapply(fits, function(x) x@opt$convergence)
out$CondNum <- sapply(fits, function(x) cn(x))
out$negLogLike <- sapply(fits, function(x) x@negLogLike)
out$nPars <- sapply(fits, function(x) length(coef(x)))
out$n <- sapply(fits, function(x) sampleSize(x))
out$AIC <- sapply(fits, function(x) x@AIC)
out$delta <- out$AIC - min(out$AIC)
out$AICwt <- exp(-out$delta / 2)
out$AICwt <- out$AICwt / sum(out$AICwt)
out$Rsq <- NA
if(!is.null(nullmod)) {
if (is.na(match(nullmod, names(fits)))) {
stop(paste("No fit named", nullmod, "was found in fits."))
}
nullmod <- fits[[nullmod]]
out$Rsq <- sapply(fits, nagR2, nullmod)
}
out <- out[order(out$AIC),]
out$cumltvWt <- cumsum(out$AICwt)
msout <- new("unmarkedModSel", Full = out,
Names = rbind(Coefs = eNames, SEs = seNames))
return(msout)
})
setAs("unmarkedModSel", "data.frame", function(from) {
out <- from@Full
out
})
setMethod("show", "unmarkedModSel", function(object)
{
out <- as(object, "data.frame")
rownames(out) <- out$model
out <- out[,c('n', 'nPars', 'AIC', 'delta', 'AICwt', 'cumltvWt', 'Rsq')]
if (all(is.na(out$Rsq))) out$Rsq <- NULL
print(format(out, digits=2, nsmall=2))
})
setMethod("coef", "unmarkedModSel", function(object)
{
coefNames <- object@Names["Coefs",]
msdf <- as(object, "data.frame")
rownames(msdf) <- msdf$model
out <- msdf[,coefNames]
out
})
setMethod("SE", "unmarkedModSel", function(obj)
{
seNames <- obj@Names["SEs",]
msdf <- as(obj, "data.frame")
rownames(msdf) <- msdf$model
out <- msdf[,seNames]
out
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.