diff --git a/DESCRIPTION b/DESCRIPTION index a251b160..31d9d7cf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.4.1.9016 -Date: 2024-08-21 +Version: 1.4.1.9017 +Date: 2024-08-22 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( @@ -47,12 +47,13 @@ LazyData: yes Collate: 'classes.R' 'unmarkedEstimate.R' 'unmarkedFrame.R' 'unmarkedFit.R' 'utils.R' 'getDesign.R' 'fitted.R' 'residuals.R' 'update.R' 'simulate.R' 'nonparboot.R' 'getP.R' + 'ranef.R' 'posteriorSamples.R' 'colext.R' 'distsamp.R' 'multinomPois.R' 'occu.R' 'occuRN.R' 'occuMulti.R' 'pcount.R' 'gmultmix.R' 'pcountOpen.R' 'gdistsamp.R' 'unmarkedFitList.R' 'unmarkedLinComb.R' - 'ranef.R' 'parboot.R' 'occuFP.R' 'gpcount.R' 'occuPEN.R' 'pcount.spHDS.R' + 'parboot.R' 'occuFP.R' 'gpcount.R' 'occuPEN.R' 'pcount.spHDS.R' 'occuMS.R' 'occuTTD.R' 'distsampOpen.R' 'multmixOpen.R' - 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'posteriorSamples.R' + 'unmarkedCrossVal.R' 'piFun.R' 'vif.R' 'makePiFun.R' 'nmixTTD.R' 'gdistremoval.R' 'plotEffects.R' diff --git a/R/IDS.R b/R/IDS.R index a8c69e6b..b0adf674 100644 --- a/R/IDS.R +++ b/R/IDS.R @@ -596,8 +596,8 @@ setMethod("nonparboot_internal", "unmarkedFitIDS", stop("Not currently supported for unmarkedFitIDS", call.=FALSE) }) -setMethod("ranef", "unmarkedFitIDS", - function(object, B = 0, keepOldSamples = TRUE, ...) +setMethod("ranef_internal", "unmarkedFitIDS", + function(object, ...) { stop("Not currently supported for unmarkedFitIDS", call.=FALSE) }) diff --git a/R/gdistremoval.R b/R/gdistremoval.R index 9c23f90f..a1fb57fa 100644 --- a/R/gdistremoval.R +++ b/R/gdistremoval.R @@ -595,7 +595,7 @@ setMethod("residuals_internal", "unmarkedFitGDR", function(object){ # ranef -setMethod("ranef", "unmarkedFitGDR", function(object){ +setMethod("ranef_internal", "unmarkedFitGDR", function(object, ...){ M <- numSites(object@data) T <- object@data@numPrimary diff --git a/R/goccu.R b/R/goccu.R index 0f80a181..685aacc6 100644 --- a/R/goccu.R +++ b/R/goccu.R @@ -208,7 +208,7 @@ setMethod("fitted_internal", "unmarkedFitGOccu", function(object){ # based on ranef for GPC -setMethod("ranef", "unmarkedFitGOccu", function(object, ...){ +setMethod("ranef_internal", "unmarkedFitGOccu", function(object, ...){ M <- numSites(object@data) JT <- obsNum(object@data) diff --git a/R/occuCOP.R b/R/occuCOP.R index 23a8f55a..70da53bf 100644 --- a/R/occuCOP.R +++ b/R/occuCOP.R @@ -458,7 +458,7 @@ setMethod("nonparboot_internal", "unmarkedFitOccuCOP", ## ranef ---- -setMethod("ranef", "unmarkedFitOccuCOP", function(object, ...) { +setMethod("ranef_internal", "unmarkedFitOccuCOP", function(object, ...) { # Sites removed (srm) and sites kept (sk) srm <- object@sitesRemoved if (length(srm) > 0) { diff --git a/R/posteriorSamples.R b/R/posteriorSamples.R index 9389e141..2d59971e 100644 --- a/R/posteriorSamples.R +++ b/R/posteriorSamples.R @@ -19,6 +19,7 @@ setMethod("posteriorSamples", "unmarkedRanef", function(object, nsims=100, ...) out <- array(NA, c(N, T, nsims)) for (n in 1:N){ + if(any(is.na(object@post[n,,]))) next for (t in 1:T){ out[n, t, ] <- sample(0:(K-1), nsims, replace=TRUE, prob=object@post[n,,t]) diff --git a/R/ranef.R b/R/ranef.R index e2f16f2a..51858e76 100644 --- a/R/ranef.R +++ b/R/ranef.R @@ -1,289 +1,270 @@ - - - - # ----------------- Empirical Bayes Methods ------------------------------ - - setGeneric("ranef", function(object, ...) standardGeneric("ranef")) - setClass("unmarkedRanef", representation(post = "array")) +# Overall exported function +setMethod("ranef", "unmarkedFit", function(object, ...){ + ranef_internal(object, ...) +}) -setClassUnion("unmarkedFitGMMorGDS", - c("unmarkedFitGMM", "unmarkedFitGDS")) +# Internal fit-type-specific function +setGeneric("ranef_internal", function(object, ...) standardGeneric("ranef_internal")) -setMethod("ranef", "unmarkedFitPCount", - function(object, ...) -{ - lam <- predict(object, type="state")[,1] # Too slow - R <- length(lam) - p <- getP(object) - K <- object@K - N <- 0:K - y <- getY(getData(object)) - srm <- object@sitesRemoved - if(length(srm) > 0) - y <- y[-object@sitesRemoved,] - post <- array(NA_real_, c(R, length(N), 1)) - colnames(post) <- N - mix <- object@mixture - for(i in 1:R) { - switch(mix, - P = f <- dpois(N, lam[i], log=TRUE), - NB = { - alpha <- exp(coef(object, type="alpha")) - f <- dnbinom(N, mu=lam[i], size=alpha, log=TRUE) - }, - ZIP = { - psi <- plogis(coef(object, type="psi")) - f <- (1-psi)*dpois(N, lam[i]) - f[1] <- psi + (1-psi)*exp(-lam[i]) - f <- log(f) - }) - g <- rep(0, K+1) - for(j in 1:ncol(y)) { - if(is.na(y[i,j]) | is.na(p[i,j])) - next - g <- g + dbinom(y[i,j], N, p[i,j], log=TRUE) - } - fudge <- exp(f+g) - post[i,,1] <- fudge / sum(fudge) - } - new("unmarkedRanef", post=post) -}) +setMethod("ranef_internal", "unmarkedFitColExt", function(object){ + data <- object@data + M <- numSites(data) + nY <- data@numPrimary + J <- obsNum(data)/nY + y <- data@y + y[y>1] <- 1 + ya <- array(y, c(M, J, nY)) + psiP <- predict(object, type="psi", level=NULL, na.rm=FALSE)$Predicted + detP <- predict(object, type="det", level=NULL, na.rm=FALSE)$Predicted + colP <- predict(object, type="col", level=NULL, na.rm=FALSE)$Predicted + extP <- predict(object, type="ext", level=NULL, na.rm=FALSE)$Predicted + detP <- array(detP, c(J, nY, M)) + colP <- matrix(colP, M, nY, byrow = TRUE) + extP <- matrix(extP, M, nY, byrow = TRUE) + ## create transition matrices (phi^T) + phis <- array(NA,c(2,2,nY-1,M)) #array of phis for each + for(i in 1:M) { + for(t in 1:(nY-1)) { + phis[,,t,i] <- matrix(c(1-colP[i,t], colP[i,t], extP[i,t], + 1-extP[i,t])) + } + } + ## first compute latent probs + x <- array(NA, c(2, nY, M)) + x[1,1,] <- 1-psiP + x[2,1,] <- psiP + for(i in 1:M) { + for(t in 2:nY) { + x[,t,i] <- (phis[,,t-1,i] %*% x[,t-1,i]) + } + } -setMethod("ranef", "unmarkedFitOccu", - function(object, ...) -{ - psi <- predict(object, type="state")[,1] - R <- length(psi) - p <- getP(object) z <- 0:1 - y <- getY(getData(object)) - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - p <- p[-object@sitesRemoved,] # temporary workaround - } - y[y>1] <- 1 - post <- array(0, c(R,2,1)) + post <- array(NA_real_, c(M, 2, nY)) colnames(post) <- z - for(i in 1:R) { - f <- dbinom(z, 1, psi[i]) - g <- rep(1, 2) - for(j in 1:ncol(y)) { - if(is.na(y[i,j]) | is.na(p[i,j])) - next - g <- g * dbinom(y[i,j], 1, z*p[i,j]) + + for(i in 1:M) { + for(t in 1:nY) { + g <- rep(1, 2) + for(j in 1:J) { + if(is.na(ya[i,j,t]) | is.na(detP[j,t,i])) + next + g <- g * dbinom(ya[i,j,t], 1, z*detP[j,t,i]) + } + tmp <- x[,t,i] * g + post[i,,t] <- tmp/sum(tmp) } - fudge <- f*g - post[i,,1] <- fudge / sum(fudge) } + new("unmarkedRanef", post=post) }) -setMethod("ranef", "unmarkedFitOccuMS", function(object, ...) -{ - - N <- numSites(object@data) - S <- object@data@numStates - - psi <- predict(object, "psi", level=NULL) - psi <- sapply(psi, function(x) x$Predicted) - z <- 0:(S-1) +# DSO and MMO +setMethod("ranef_internal", "unmarkedFitDailMadsen", function(object, ...){ + dyn <- object@dynamics + formlist <- object@formlist + formula <- as.formula(paste(unlist(formlist), collapse=" ")) + D <- getDesign(object@data, formula, na.rm=FALSE) + delta <- D$delta + deltamax <- max(delta, na.rm=TRUE) + if(!.hasSlot(object, "immigration")){ #For backwards compatibility + imm <- FALSE + } else { + imm <- object@immigration + } - p_all <- getP(object) - y <- getY(getData(object)) + #TODO: adjust if output = "density" + lam <- predict(object, type="lambda",level=NULL, na.rm=FALSE)$Predicted # Slow, use D$Xlam instead - post <- array(0, c(N,S,1)) - colnames(post) <- z + R <- length(lam) + T <- object@data@numPrimary - if(object@parameterization == "multinomial"){ + p <- getP(object) - psi <- cbind(1-rowSums(psi), psi) + K <- object@K + N <- 0:K + y <- getY(getData(object)) + J <- ncol(y)/T + if(dyn != "notrend") { + gam <- predict(object, type="gamma", level=NULL, na.rm=FALSE)$Predicted + gam <- matrix(gam, R, T-1, byrow=TRUE) + } + if(!identical(dyn, "trend")) { + om <- predict(object, type="omega", level=NULL, na.rm=FALSE)$Predicted + om <- matrix(om, R, T-1, byrow=TRUE) + } else { + om <- matrix(0, R, T-1) + } + if(imm) { + iota <- predict(object, type="iota",level=NULL,na.rm=FALSE)$Predicted + iota <- matrix(iota, R, T-1, byrow=TRUE) + } else { + iota <- matrix(0, R, T-1) + } + ya <- array(y, c(R, J, T)) + pa <- array(p, c(R, J, T)) + post <- array(NA_real_, c(R, length(N), T)) + colnames(post) <- N + mix <- object@mixture + if(dyn=="notrend") + gam <- lam*(1-om) - guide <- matrix(NA,nrow=S,ncol=S) - guide <- lower.tri(guide,diag=TRUE) - guide[,1] <- FALSE - guide <- which(guide,arr.ind=TRUE) - for (i in 1:N){ - f <- psi[i,] - g <- rep(1, S) - p_raw <- sapply(p_all, function(x) x[i,]) - for (j in 1:nrow(p_raw)){ - if(any(is.na(p_raw[j,])) | is.na(y[i,j])) next - sdp <- matrix(0, nrow=S, ncol=S) - sdp[guide] <- p_raw[j,] - sdp[,1] <- 1 - rowSums(sdp) - for (s in 1:S){ - g[s] <- g[s] * sdp[s, (y[i,j]+1)] + if(dyn %in% c("constant", "notrend")) { + tp <- function(N0, N1, gam, om, iota) { + c <- 0:min(N0, N1) + sum(dbinom(c, N0, om) * dpois(N1-c, gam)) + } + } else if(dyn=="autoreg") { + tp <- function(N0, N1, gam, om, iota) { + c <- 0:min(N0, N1) + sum(dbinom(c, N0, om) * dpois(N1-c, gam*N0 + iota)) + } + } else if(dyn=="trend") { + tp <- function(N0, N1, gam, om, iota) { + dpois(N1, gam*N0 + iota) + } + } else if(dyn=="ricker") { + tp <- function(N0, N1, gam, om, iota) { + dpois(N1, N0*exp(gam*(1-N0/om)) + iota) + } + } else if(dyn=="gompertz") { + tp <- function(N0, N1, gam, om, iota) { + dpois(N1, N0*exp(gam*(1-log(N0 + 1)/log(om + 1))) + iota) } - } - fudge <- f*g - post[i,,1] <- fudge / sum(fudge) } + for(i in 1:R) { + P <- matrix(1, K+1, K+1) + switch(mix, + P = g2 <- dpois(N, lam[i]), + NB = { + alpha <- exp(coef(object, type="alpha")) + g2 <- dnbinom(N, mu=lam[i], size=alpha) + }, + ZIP = { + psi <- plogis(coef(object, type="psi")) + g2 <- (1-psi)*dpois(N, lam[i]) + g2[1] <- psi + (1-psi)*exp(-lam[i]) + }) - } else if(object@parameterization == "condbinom"){ - - psi <- cbind(1-psi[,1], psi[,1]*(1-psi[,2]), psi[,1]*psi[,2]) + #DETECTION MODEL + g1 <- rep(0, K+1) + cp <- pa[i,,1] + cp_na <- is.na(cp) + ysub <- ya[i,,1] + ysub[cp_na] <- NA + cp <- c(cp, 1-sum(cp, na.rm=TRUE)) + sumy <- sum(ysub, na.rm=TRUE) - for (i in 1:N){ - f <- psi[i,] - g <- rep(1, S) - p_raw <- sapply(p_all, function(x) x[i,]) - for (j in 1:nrow(p_raw)){ - probs <- p_raw[j,] - if(any(is.na(probs)) | is.na(y[i,j])) next - sdp <- matrix(0, nrow=S, ncol=S) - sdp[1,1] <- 1 - sdp[2,1:2] <- c(1-probs[1], probs[1]) - sdp[3,] <- c(1-probs[2], probs[2]*(1-probs[3]), probs[2]*probs[3]) - for (s in 1:S){ - g[s] <- g[s] * sdp[s, (y[i,j]+1)] - } - } - fudge <- f*g - post[i,,1] <- fudge / sum(fudge) - } - } + is_na <- c(is.na(ysub), FALSE) | is.na(cp) - new("unmarkedRanef", post=post) + if(all(is.na(ysub))){ + post[i,,1] <- NA + } else { -}) + for(k in sumy:K){ + yit <- c(ysub, k-sumy) + g1[k+1] <- dmultinom(yit[!is_na], k, cp[!is_na]) + } -setMethod("ranef", "unmarkedFitOccuFP", function(object, na.rm = FALSE) -{ - cat("ranef is not implemented for occuFP at this time") -}) + g1g2 <- g1*g2 + post[i,,1] <- g1g2 / sum(g1g2) + } -setMethod("ranef", "unmarkedFitOccuMulti", function(object, species, na.rm = FALSE) -{ - if(missing(species)){ - stop("Must specify species name or index in arguments (e.g. species = 1)") - } - species <- name_to_ind(species, names(object@data@ylist)) + for(t in 2:T) { + if(!is.na(gam[i,t-1]) & !is.na(om[i,t-1])) { + for(n0 in N) { + for(n1 in N) { + P[n0+1, n1+1] <- tp(n0, n1, gam[i,t-1], om[i,t-1], iota[i,t-1]) + } + } + } + delta.it <- delta[i,t-1] + if(delta.it > 1) { + P1 <- P + for(d in 2:delta.it) { + P <- P %*% P1 + } + } - psi <- predict(object, type="state", level=NULL, species=species)$Predicted - R <- length(psi) - p <- getP(object)[[species]] - z <- 0:1 - y <- object@data@ylist[[species]] - post <- array(0, c(R,2,1)) - colnames(post) <- z - for(i in 1:R) { - f <- dbinom(z, 1, psi[i]) - g <- rep(1, 2) - for(j in 1:ncol(y)) { - if(is.na(y[i,j]) | is.na(p[i,j])) - next - g <- g * dbinom(y[i,j], 1, z*p[i,j]) - } - fudge <- f*g - post[i,,1] <- fudge / sum(fudge) - } - new("unmarkedRanef", post=post) -}) + #DETECTION MODEL + g1 <- rep(0, K+1) + cp <- pa[i,,t] + cp_na <- is.na(cp) + ysub <- ya[i,,t] + ysub[cp_na] <- NA + cp <- c(cp, 1-sum(cp, na.rm=TRUE)) + sumy <- sum(ysub, na.rm=TRUE) + is_na <- c(is.na(ysub), FALSE) | is.na(cp) + if(all(is.na(ysub))){ + post[i,,t] <- NA + } else { + for(k in sumy:K){ + yit <- c(ysub, k-sumy) + g1[k+1] <- dmultinom(yit[!is_na], k, cp[!is_na]) + } -setMethod("ranef", "unmarkedFitOccuRN", - function(object, K, ...) -{ - if(missing(K)) { - warning("You did not specify K, the maximum value of N, so it was set to 50") - K <- 50 - } - lam <- predict(object, type="state")[,1] # Too slow - R <- length(lam) - r <- getP(object) - N <- 0:K - y <- getY(getData(object)) - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - r <- r[-object@sitesRemoved,] # temporary workaround - } - y[y>1] <- 1 - post <- array(NA_real_, c(R, length(N), 1)) - colnames(post) <- N - for(i in 1:R) { - f <- dpois(N, lam[i]) - g <- rep(1, K+1) - for(j in 1:ncol(y)) { - if(is.na(y[i,j]) | is.na(r[i,j])) - next - p.ijn <- 1 - (1-r[i,j])^N - g <- g * dbinom(y[i,j], 1, p.ijn) - } - fudge <- f*g - post[i,,1] <- fudge / sum(fudge) + g <- colSums(P * post[i,,t-1]) * g1 + post[i,,t] <- g / sum(g) + } + } } new("unmarkedRanef", post=post) }) - - - -setMethod("ranef", "unmarkedFitMPois", - function(object, K, ...) -{ +setMethod("ranef_internal", "unmarkedFitDS", function(object, ...){ y <- getY(getData(object)) cp <- getP(object) - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - cp <- cp[-object@sitesRemoved,] # temporary workaround - } - if(missing(K)) { - warning("You did not specify K, the maximum value of N, so it was set to max(y)+50") - K <- max(y, na.rm=TRUE)+50 + K <- list(...)$K + if(is.null(K)) { + warning("You did not specify K, the maximum value of N, so it was set to max(y)+50") + K <- max(y, na.rm=TRUE)+50 } - - lam <- predict(object, type="state")[,1] + lam <- predict(object, type="state", level=NULL, na.rm=FALSE)$Predicted R <- length(lam) - cp <- cbind(cp, 1-rowSums(cp, na.rm=TRUE)) + J <- ncol(y) + if(identical(object@output, "density")) { + A <- get_ds_area(object@data, object@unitsOut) + lam <- lam*A + } + cp <- cbind(cp, 1-rowSums(cp)) N <- 0:K - post <- array(NA, c(R, K+1, 1)) + post <- array(0, c(R, K+1, 1)) colnames(post) <- N for(i in 1:R) { f <- dpois(N, lam[i]) g <- rep(1, K+1) - yi <- y[i,] - cpi <- cp[i,] if(any(is.na(y[i,])) | any(is.na(cp[i,]))){ - # This only handles cases when all NAs are at the end of the count vector - # Not when they are interspersed. I don't know how to handle that situation - if(object@data@piFun == "removalPiFun"){ - warning("NAs in counts and/or covariates for site ",i, ". Keeping only counts before the first NA", call.=FALSE) - notNA <- min(which(is.na(y[i,]))[1], which(is.na(cp[i,]))[1], na.rm=TRUE) - 1 - yi <- yi[c(1:notNA)] - cpi <- cpi[c(1:notNA, length(cpi))] - } else { + post[i,,1] <- NA next - } } for(k in 1:(K+1)) { + yi <- y[i,] ydot <- N[k] - sum(yi) if(ydot<0) { g[k] <- 0 next } - yik <- c(yi, ydot) - g[k] <- g[k] * dmultinom(yik, size=N[k], prob=cpi) + yi <- c(yi, ydot) + g[k] <- g[k] * dmultinom(yi, size=N[k], prob=cp[i,]) } fudge <- f*g post[i,,1] <- fudge / sum(fudge) @@ -292,70 +273,74 @@ setMethod("ranef", "unmarkedFitMPois", }) +setMethod("ranef_internal", "unmarkedFitOccu", function(object, ...){ + psi <- predict(object, type="state", level=NULL, na.rm=FALSE)$Predicted + R <- length(psi) + p <- getP(object) + z <- 0:1 + y <- getY(getData(object)) + y[y>1] <- 1 + post <- array(NA, c(R,2,1)) + colnames(post) <- z + for(i in 1:R) { + if(all(is.na(y[i,]))) next + f <- dbinom(z, 1, psi[i]) + g <- rep(1, 2) + for(j in 1:ncol(y)) { + if(is.na(y[i,j]) | is.na(p[i,j])) + next + g <- g * dbinom(y[i,j], 1, z*p[i,j]) + } + fudge <- f*g + post[i,,1] <- fudge / sum(fudge) + } + new("unmarkedRanef", post=post) +}) - - - -setMethod("ranef", "unmarkedFitDS", - function(object, K, ...) -{ +setMethod("ranef_internal", "unmarkedFitMPois", function(object, ...){ y <- getY(getData(object)) cp <- getP(object) - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - cp <- cp[-object@sitesRemoved,] # temporary workaround - } - if(missing(K)) { + K <- list(...)$K + if(is.null(K)) { warning("You did not specify K, the maximum value of N, so it was set to max(y)+50") K <- max(y, na.rm=TRUE)+50 } - lam <- predict(object, type="state")[,1] + + lam <- predict(object, type="state", level=NULL, na.rm=FALSE)$Predicted R <- length(lam) - J <- ncol(y) - survey <- object@data@survey - tlength <- object@data@tlength - db <- object@data@dist.breaks - w <- diff(db) - unitsIn <- object@data@unitsIn - unitsOut <- object@unitsOut - if(identical(object@output, "density")) { - a <- matrix(NA, R, J) - switch(survey, line = { - for (i in 1:R) { - a[i, ] <- tlength[i] * w - } - }, point = { - for (i in 1:R) { - a[i, 1] <- pi * db[2]^2 - for (j in 2:J) - a[i, j] <- pi * db[j + 1]^2 - sum(a[i, 1:(j - 1)]) - } - }) - switch(survey, line = A <- rowSums(a) * 2, point = A <- rowSums(a)) - switch(unitsIn, m = A <- A/1e+06, km = A <- A) - switch(unitsOut, ha = A <- A * 100, kmsq = A <- A) - lam <- lam*A - } - cp <- cbind(cp, 1-rowSums(cp)) + all_na <- apply(cp, 1, function(x) all(is.na(x))) + cp <- cbind(cp, 1-rowSums(cp, na.rm=TRUE)) + cp[all_na,ncol(cp)] <- NA N <- 0:K - post <- array(0, c(R, K+1, 1)) + post <- array(NA, c(R, K+1, 1)) colnames(post) <- N for(i in 1:R) { f <- dpois(N, lam[i]) g <- rep(1, K+1) - if(any(is.na(y[i,])) | any(is.na(cp[i,]))) + yi <- y[i,] + cpi <- cp[i,] + if(any(is.na(y[i,])) | any(is.na(cp[i,]))){ + # This only handles cases when all NAs are at the end of the count vector + # Not when they are interspersed. I don't know how to handle that situation + if(object@data@piFun == "removalPiFun"){ + if(all(is.na(cp[i,]))) next + warning("NAs in counts and/or covariates for site ",i, ". Keeping only counts before the first NA", call.=FALSE) + notNA <- min(which(is.na(y[i,]))[1], which(is.na(cp[i,]))[1], na.rm=TRUE) - 1 + yi <- yi[c(1:notNA)] + cpi <- cpi[c(1:notNA, length(cpi))] + } else { next + } + } for(k in 1:(K+1)) { - yi <- y[i,] ydot <- N[k] - sum(yi) if(ydot<0) { g[k] <- 0 next } - yi <- c(yi, ydot) - g[k] <- g[k] * dmultinom(yi, size=N[k], prob=cp[i,]) + yik <- c(yi, ydot) + g[k] <- g[k] * dmultinom(yik, size=N[k], prob=cpi) } fudge <- f*g post[i,,1] <- fudge / sum(fudge) @@ -364,89 +349,40 @@ setMethod("ranef", "unmarkedFitDS", }) +setMethod("ranef_internal", "unmarkedFitOccuFP", function(object, ...){ + stop("ranef is not implemented for occuFP at this time", call.=FALSE) +}) - - - -setMethod("ranef", "unmarkedFitGMMorGDS", - function(object, ...) -{ - +# Function that works for both GMM and GDS +# Avoiding class union that doesn't work for some reason +ranef_GMM_GDS <- function(object, ...){ data <- object@data - D <- getDesign(data, object@formula) - - Xlam <- D$Xlam - Xphi <- D$Xphi - Xdet <- D$Xdet - y <- D$y # MxJT - Xlam.offset <- D$Xlam.offset - Xphi.offset <- D$Xphi.offset - Xdet.offset <- D$Xdet.offset - if(is.null(Xlam.offset)) Xlam.offset <- rep(0, nrow(Xlam)) - if(is.null(Xphi.offset)) Xphi.offset <- rep(0, nrow(Xphi)) - if(is.null(Xdet.offset)) Xdet.offset <- rep(0, nrow(Xdet)) - beta.lam <- coef(object, type="lambda") - beta.phi <- coef(object, type="phi") - beta.det <- coef(object, type="det") - - lambda <- exp(Xlam %*% beta.lam + Xlam.offset) - if(is.null(beta.phi)) - phi <- rep(1, nrow(Xphi)) - else - phi <- plogis(Xphi %*% beta.phi + Xphi.offset) + y <- getY(data) + nSites <- numSites(data) + T <- data@numPrimary + R <- numY(data) / T + J <- obsNum(data) / T - cp <- getP(object) - srm <- object@sitesRemoved - if(length(srm) > 0){ - cp <- cp[-object@sitesRemoved,] # temporary workaround + lambda <- predict(object, type="lambda", level=NULL, na.rm=FALSE)$Predicted + if(T == 1){ + phi <- rep(1, nSites) + } else { + phi <- predict(object, type="phi", level=NULL, na.rm=FALSE)$Predicted } + + cp <- getP(object) cp[is.na(y)] <- NA K <- object@K M <- 0:K - nSites <- nrow(y) - T <- data@numPrimary - R <- numY(data) / T - J <- obsNum(data) / T - phi <- matrix(phi, nSites, byrow=TRUE) if(inherits(object, "unmarkedFitGDS")) { if(identical(object@output, "density")) { - survey <- object@data@survey - tlength <- object@data@tlength - unitsIn <- object@data@unitsIn - unitsOut <- object@unitsOut - db <- object@data@dist.breaks - w <- diff(db) - u <- a <- matrix(NA, nSites, R) - switch(survey, - line = { - for(i in 1:nSites) { - a[i,] <- tlength[i] * w - u[i,] <- a[i,] / sum(a[i,]) - } - }, - point = { - for(i in 1:nSites) { - a[i, 1] <- pi*db[2]^2 - for(j in 2:R) - a[i, j] <- pi*db[j+1]^2 - sum(a[i, 1:(j-1)]) - u[i,] <- a[i,] / sum(a[i,]) - } - }) - switch(survey, - line = A <- rowSums(a) * 2, - point = A <- rowSums(a)) - switch(unitsIn, - m = A <- A / 1e6, - km = A <- A) - switch(unitsOut, - ha = A <- A * 100, - kmsq = A <- A) - lambda <- lambda*A # Density to abundance + A <- get_ds_area(data, object@unitsOut) + lambda <- lambda*A # Density to abundance } } @@ -454,7 +390,7 @@ setMethod("ranef", "unmarkedFitGMMorGDS", ya <- array(y, c(nSites,R,T)) # ym <- apply(ya, c(1,3), sum, na.rm=TRUE) - post <- array(0, c(nSites, K+1, 1)) + post <- array(NA_real_, c(nSites, K+1, 1)) colnames(post) <- M mix <- object@mixture if(identical(mix, "NB")){ @@ -463,13 +399,14 @@ setMethod("ranef", "unmarkedFitGMMorGDS", psi <- plogis(coef(object, type="psi")) } for(i in 1:nSites) { + if(all(is.na(ya[i,,]))) next switch(mix, P = f <- dpois(M, lambda[i]), ZIP = f <- dzip(M, lambda[i], psi), NB = f <- dnbinom(M, mu=lambda[i], size=alpha)) g <- rep(1, K+1) # outside t loop for(t in 1:T) { - if(all(is.na(ya[i,,t])) | is.na(phi[i,t])) + if(any(is.na(ya[i,,t])) | any(is.na(cpa[i,,t])) | is.na(phi[i,t])) next for(k in 1:(K+1)) { y.it <- ya[i,,t] @@ -490,50 +427,33 @@ setMethod("ranef", "unmarkedFitGMMorGDS", post[i,,1] <- fudge/sum(fudge) } new("unmarkedRanef", post=post) -}) - +} +setMethod("ranef_internal", "unmarkedFitGDS", function(object, ...){ + ranef_GMM_GDS(object, ...) +}) +setMethod("ranef_internal", "unmarkedFitGMM", function(object, ...){ + ranef_GMM_GDS(object, ...) +}) -setMethod("ranef", "unmarkedFitGPC", - function(object, ...) -{ +setMethod("ranef_internal", "unmarkedFitGPC", function(object, ...){ data <- object@data - D <- getDesign(data, object@formula) - - Xlam <- D$Xlam - Xphi <- D$Xphi - Xdet <- D$Xdet - y <- D$y # MxJT - Xlam.offset <- D$Xlam.offset - Xphi.offset <- D$Xphi.offset - Xdet.offset <- D$Xdet.offset - if(is.null(Xlam.offset)) Xlam.offset <- rep(0, nrow(Xlam)) - if(is.null(Xphi.offset)) Xphi.offset <- rep(0, nrow(Xphi)) - if(is.null(Xdet.offset)) Xdet.offset <- rep(0, nrow(Xdet)) - - beta.lam <- coef(object, type="lambda") - beta.phi <- coef(object, type="phi") - beta.det <- coef(object, type="det") - - R <- nrow(y) + R <- numSites(data) T <- data@numPrimary + y <- getY(data) J <- ncol(y) / T - lambda <- exp(Xlam %*% beta.lam + Xlam.offset) - if(is.null(beta.phi)) - phi <- rep(1, nrow(Xphi)) + lambda <- predict(object, type="lambda", level=NULL, na.rm=FALSE)$Predicted + if(T == 1) + phi <- rep(1, R) else - phi <- plogis(Xphi %*% beta.phi + Xphi.offset) + phi <- predict(object, type="phi", level=NULL, na.rm=FALSE)$Predicted phi <- matrix(phi, R, byrow=TRUE) p <- getP(object) - srm <- object@sitesRemoved - if(length(srm) > 0){ - p <- p[-object@sitesRemoved,] # temporary workaround - } p[is.na(y)] <- NA pa <- array(p, c(R,J,T)) ya <- array(y, c(R,J,T)) @@ -542,7 +462,7 @@ setMethod("ranef", "unmarkedFitGPC", M <- N <- 0:K lM <- K+1 - post <- array(0, c(R, K+1, 1)) + post <- array(NA_real_, c(R, K+1, 1)) colnames(post) <- M mix <- object@mixture if(identical(mix, "NB")) @@ -551,6 +471,7 @@ setMethod("ranef", "unmarkedFitGPC", psi <- plogis(coef(object, type="psi")) for(i in 1:R) { + if(all(is.na(ya[i,,]))) next switch(mix, P = f <- dpois(M, lambda[i]), NB = f <- dnbinom(M, mu=lambda[i], size=alpha), @@ -588,224 +509,208 @@ setMethod("ranef", "unmarkedFitGPC", }) +setMethod("ranef_internal", "unmarkedFitNmixTTD", function(object, ...){ + M <- nrow(object@data@y) + J <- ncol(object@data@y) + tdist <- ifelse("shape" %in% names(object@estimates), "weibull", "exp") + mix <- ifelse("alpha" %in% names(object@estimates), "NB", "P") + K <- object@K + yvec <- as.numeric(t(object@data@y)) + removed <- object@sitesRemoved + naflag <- as.numeric(is.na(yvec)) + surveyLength <- object@data@surveyLength + if(length(removed>0)) surveyLength <- surveyLength[-removed,] + ymax <- as.numeric(t(surveyLength)) + delta <- as.numeric(yvec1] <- 1 - ya <- array(y, c(M, J, nY)) - - psiP <- plogis(W.i %*% psiParms) - detP <- plogis(V.itj %*% detParms) - colP <- plogis(X.it.gam %*% colParms) - extP <- plogis(X.it.eps %*% extParms) + if(mix == "P"){ + pK <- sapply(0:K, function(k) dpois(k, abun)) + } else { + alpha <- exp(coef(object, "alpha")) + pK <- sapply(0:K, function(k) dnbinom(k, mu=abun, size = alpha)) + } - detP <- array(detP, c(J, nY, M)) - colP <- matrix(colP, M, nY, byrow = TRUE) - extP <- matrix(extP, M, nY, byrow = TRUE) + if(tdist=='weibull'){ + shape <- exp(coef(object, "shape")) - ## create transition matrices (phi^T) - phis <- array(NA,c(2,2,nY-1,M)) #array of phis for each - for(i in 1:M) { - for(t in 1:(nY-1)) { - phis[,,t,i] <- matrix(c(1-colP[i,t], colP[i,t], extP[i,t], - 1-extP[i,t])) - } - } + e_lamt <- sapply(0:K, function(k){ + lamK <- k*lam + ( shape*lamK*(lamK*yvec)^(shape-1) )^delta * exp(-1*(lamK*yvec)^shape) + }) - ## first compute latent probs - x <- array(NA, c(2, nY, M)) - x[1,1,] <- 1-psiP - x[2,1,] <- psiP - for(i in 1:M) { - for(t in 2:nY) { - x[,t,i] <- (phis[,,t-1,i] %*% x[,t-1,i]) - } - } + } else { + #Exponential + e_lamt <- sapply(0:K, function(k) (lam*k)^delta * exp(-lam*k*yvec)) + } - z <- 0:1 - post <- array(NA_real_, c(M, 2, nY)) - colnames(post) <- z + post <- array(NA, c(M, K+1, 1)) + colnames(post) <- 0:K + ystart <- 1 + for (m in 1:M){ - for(i in 1:M) { - for(t in 1:nY) { - g <- rep(1, 2) - for(j in 1:J) { - if(is.na(ya[i,j,t]) | is.na(detP[j,t,i])) - next - g <- g * dbinom(ya[i,j,t], 1, z*detP[j,t,i]) - } - tmp <- x[,t,i] * g - post[i,,t] <- tmp/sum(tmp) - } + yend <- ystart+J-1 + pT <- rep(NA,length=K+1) + pT[1] <- 1 - max(delta[ystart:yend], na.rm=T) + for (k in 1:K){ + elamt_sub <- e_lamt[ystart:yend, k+1] + pT[k+1] <- prod(elamt_sub[!is.na(elamt_sub)]) } + ystart <- ystart + J + probs <- pK[m,] * pT + post[m,,1] <- probs / sum(probs) + } - new("unmarkedRanef", post=post) + new("unmarkedRanef", post=post) }) +setMethod("ranef_internal", "unmarkedFitOccuMS", function(object, ...){ + N <- numSites(object@data) + S <- object@data@numStates + psi <- predict(object, "psi", level=NULL) + psi <- sapply(psi, function(x) x$Predicted) + z <- 0:(S-1) + p_all <- getP(object) + y <- getY(getData(object)) + post <- array(NA_real_, c(N,S,1)) + colnames(post) <- z -setMethod("ranef", "unmarkedFitPCO", - function(object, ...) -{ -# browser() - dyn <- object@dynamics - formlist <- object@formlist - formula <- as.formula(paste(unlist(formlist), collapse=" ")) - D <- getDesign(object@data, formula) - delta <- D$delta - deltamax <- max(delta, na.rm=TRUE) - if(!.hasSlot(object, "immigration")) #For backwards compatibility - imm <- FALSE - else - imm <- object@immigration + if(object@parameterization == "multinomial"){ - lam <- predict(object, type="lambda")[,1] # Slow, use D$Xlam instead - R <- length(lam) - T <- object@data@numPrimary - p <- getP(object) - K <- object@K - N <- 0:K - y <- getY(getData(object)) - J <- ncol(y)/T - if(dyn != "notrend") { - gam <- predict(object, type="gamma")[,1] - gam <- matrix(gam, R, T-1, byrow=TRUE) - } - if(!identical(dyn, "trend")) { - om <- predict(object, type="omega")[,1] - om <- matrix(om, R, T-1, byrow=TRUE) + psi <- cbind(1-rowSums(psi), psi) + + guide <- matrix(NA,nrow=S,ncol=S) + guide <- lower.tri(guide,diag=TRUE) + guide[,1] <- FALSE + guide <- which(guide,arr.ind=TRUE) + for (i in 1:N){ + if(all(is.na(y[i,]))) next + f <- psi[i,] + g <- rep(1, S) + p_raw <- sapply(p_all, function(x) x[i,]) + for (j in 1:nrow(p_raw)){ + if(any(is.na(p_raw[j,])) | is.na(y[i,j])) next + sdp <- matrix(0, nrow=S, ncol=S) + sdp[guide] <- p_raw[j,] + sdp[,1] <- 1 - rowSums(sdp) + for (s in 1:S){ + g[s] <- g[s] * sdp[s, (y[i,j]+1)] + } + } + fudge <- f*g + post[i,,1] <- fudge / sum(fudge) } - else - om <- matrix(0, R, T-1) - if(imm) { - iota <- predict(object, type="iota")[,1] - iota <- matrix(iota, R, T-1, byrow=TRUE) + + } else if(object@parameterization == "condbinom"){ + + psi <- cbind(1-psi[,1], psi[,1]*(1-psi[,2]), psi[,1]*psi[,2]) + + for (i in 1:N){ + if(all(is.na(y[i,]))) next + f <- psi[i,] + g <- rep(1, S) + p_raw <- sapply(p_all, function(x) x[i,]) + for (j in 1:nrow(p_raw)){ + probs <- p_raw[j,] + if(any(is.na(probs)) | is.na(y[i,j])) next + sdp <- matrix(0, nrow=S, ncol=S) + sdp[1,1] <- 1 + sdp[2,1:2] <- c(1-probs[1], probs[1]) + sdp[3,] <- c(1-probs[2], probs[2]*(1-probs[3]), probs[2]*probs[3]) + for (s in 1:S){ + g[s] <- g[s] * sdp[s, (y[i,j]+1)] + } + } + fudge <- f*g + post[i,,1] <- fudge / sum(fudge) } - else - iota <- matrix(0, R, T-1) - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - p <- p[-object@sitesRemoved,] # temporary workaround + } + + new("unmarkedRanef", post=post) + +}) + + +setMethod("ranef_internal", "unmarkedFitOccuMulti", function(object, ...){ + + species <- list(...)$species + sp_names <- names(getData(object)@ylist) + if(!is.null(species)){ + species <- name_to_ind(species, sp_names) + } else { + species <- sp_names } - ya <- array(y, c(R, J, T)) - pa <- array(p, c(R, J, T)) - post <- array(NA_real_, c(R, length(N), T)) - colnames(post) <- N - mix <- object@mixture - if(dyn=="notrend") - gam <- lam*(1-om) - if(dyn %in% c("constant", "notrend")) { - tp <- function(N0, N1, gam, om, iota) { - c <- 0:min(N0, N1) - sum(dbinom(c, N0, om) * dpois(N1-c, gam)) - } - } else if(dyn=="autoreg") { - tp <- function(N0, N1, gam, om, iota) { - c <- 0:min(N0, N1) - sum(dbinom(c, N0, om) * dpois(N1-c, gam*N0 + iota)) - } - } else if(dyn=="trend") { - tp <- function(N0, N1, gam, om, iota) { - dpois(N1, gam*N0 + iota) - } - } else if(dyn=="ricker") { - tp <- function(N0, N1, gam, om, iota) { - dpois(N1, N0*exp(gam*(1-N0/om)) + iota) - } - } else if(dyn=="gompertz") { - tp <- function(N0, N1, gam, om, iota) { - dpois(N1, N0*exp(gam*(1-log(N0 + 1)/log(om + 1))) + iota) + out <- lapply(species, function(s){ + psi <- predict(object, type="state", level=NULL, species=s)$Predicted + R <- length(psi) + p <- getP(object)[[s]] + z <- 0:1 + y <- object@data@ylist[[s]] + post <- array(NA_real_, c(R,2,1)) + colnames(post) <- z + for(i in 1:R) { + if(all(is.na(y[i,]))) next + f <- dbinom(z, 1, psi[i]) + g <- rep(1, 2) + for(j in 1:ncol(y)) { + if(is.na(y[i,j]) | is.na(p[i,j])) + next + g <- g * dbinom(y[i,j], 1, z*p[i,j]) } + fudge <- f*g + post[i,,1] <- fudge / sum(fudge) + } + new("unmarkedRanef", post=post) + }) + names(out) <- species + if(length(out) == 1) out <- out[[1]] + out +}) + + +setMethod("ranef_internal", "unmarkedFitOccuRN", function(object, ...){ + K <- list(...)$K + if(is.null(K)) { + warning("You did not specify K, the maximum value of N, so it was set to 50") + K <- 50 } + lam <- predict(object, type="state", level=NULL, na.rm=FALSE)$Predicted # Too slow + R <- length(lam) + r <- getP(object) + N <- 0:K + y <- getY(getData(object)) + y[y>1] <- 1 + post <- array(NA_real_, c(R, length(N), 1)) + colnames(post) <- N for(i in 1:R) { - P <- matrix(1, K+1, K+1) - switch(mix, - P = g2 <- dpois(N, lam[i]), - NB = { - alpha <- exp(coef(object, type="alpha")) - g2 <- dnbinom(N, mu=lam[i], size=alpha) - }, - ZIP = { - psi <- plogis(coef(object, type="psi")) - g2 <- (1-psi)*dpois(N, lam[i]) - g2[1] <- psi + (1-psi)*exp(-lam[i]) - }) - g1 <- rep(1, K+1) - for(j in 1:J) { - if(is.na(ya[i,j,1]) | is.na(pa[i,j,1])) + if(all(is.na(y[i,]))) next + f <- dpois(N, lam[i]) + g <- rep(1, K+1) + for(j in 1:ncol(y)) { + if(is.na(y[i,j]) | is.na(r[i,j])) next - g1 <- g1 * dbinom(ya[i,j,1], N, pa[i,j,1]) - } - g1g2 <- g1*g2 - post[i,,1] <- g1g2 / sum(g1g2) - for(t in 2:T) { - if(!is.na(gam[i,t-1]) & !is.na(om[i,t-1])) { - for(n0 in N) { - for(n1 in N) { - P[n0+1, n1+1] <- tp(n0, n1, gam[i,t-1], om[i,t-1], iota[i,t-1]) - } - } - } - delta.it <- delta[i,t-1] - if(delta.it > 1) { - P1 <- P - for(d in 2:delta.it) { - P <- P %*% P1 - } - } - g1 <- rep(1, K+1) - for(j in 1:J) { - if(is.na(ya[i,j,t]) | is.na(pa[i,j,t])) - next - g1 <- g1 * dbinom(ya[i,j,t], N, pa[i,j,t]) - } - g <- colSums(P * post[i,,t-1]) * g1 - post[i,,t] <- g / sum(g) + p.ijn <- 1 - (1-r[i,j])^N + g <- g * dbinom(y[i,j], 1, p.ijn) } + fudge <- f*g + post[i,,1] <- fudge / sum(fudge) } new("unmarkedRanef", post=post) }) -setMethod("ranef", "unmarkedFitOccuTTD", - function(object, ...) -{ - +setMethod("ranef_internal", "unmarkedFitOccuTTD", function(object, ...){ N <- nrow(object@data@y) T <- object@data@numPrimary J <- ncol(object@data@y)/T @@ -866,55 +771,81 @@ setMethod("ranef", "unmarkedFitOccuTTD", }) -# DSO and MMO -setMethod("ranef", "unmarkedFitDailMadsen", - function(object, ...) -{ +setMethod("ranef_internal", "unmarkedFitPCount", function(object, ...){ + lam <- predict(object, type="state", level=NULL, na.rm=FALSE)$Predicted + R <- length(lam) + p <- getP(object) + K <- object@K + N <- 0:K + y <- getY(getData(object)) + post <- array(NA_real_, c(R, length(N), 1)) + colnames(post) <- N + mix <- object@mixture + for(i in 1:R) { + if(all(is.na(y[i,]))) next + switch(mix, + P = f <- dpois(N, lam[i], log=TRUE), + NB = { + alpha <- exp(coef(object, type="alpha")) + f <- dnbinom(N, mu=lam[i], size=alpha, log=TRUE) + }, + ZIP = { + psi <- plogis(coef(object, type="psi")) + f <- (1-psi)*dpois(N, lam[i]) + f[1] <- psi + (1-psi)*exp(-lam[i]) + f <- log(f) + }) + g <- rep(0, K+1) + for(j in 1:ncol(y)) { + if(is.na(y[i,j]) | is.na(p[i,j])) + next + g <- g + dbinom(y[i,j], N, p[i,j], log=TRUE) + } + fudge <- exp(f+g) + post[i,,1] <- fudge / sum(fudge) + } + new("unmarkedRanef", post=post) +}) + + +setMethod("ranef_internal", "unmarkedFitPCO", function(object, ...){ dyn <- object@dynamics + formlist <- object@formlist formula <- as.formula(paste(unlist(formlist), collapse=" ")) - D <- getDesign(object@data, formula) + D <- getDesign(object@data, formula, na.rm=FALSE) delta <- D$delta deltamax <- max(delta, na.rm=TRUE) - if(!.hasSlot(object, "immigration")){ #For backwards compatibility + if(!.hasSlot(object, "immigration")) #For backwards compatibility imm <- FALSE - } else { + else imm <- object@immigration - } - - #TODO: adjust if output = "density" - lam <- predict(object, type="lambda")[,1] # Slow, use D$Xlam instead + lam <- predict(object, type="lambda", level=NULL, na.rm=FALSE)$Predicted # Slow, use D$Xlam instead R <- length(lam) T <- object@data@numPrimary - p <- getP(object) - K <- object@K N <- 0:K y <- getY(getData(object)) J <- ncol(y)/T if(dyn != "notrend") { - gam <- predict(object, type="gamma")[,1] + gam <- predict(object, type="gamma", level=NULL, na.rm=FALSE)$Predicted gam <- matrix(gam, R, T-1, byrow=TRUE) } if(!identical(dyn, "trend")) { - om <- predict(object, type="omega")[,1] + om <- predict(object, type="omega", level=NULL, na.rm=FALSE)$Predicted om <- matrix(om, R, T-1, byrow=TRUE) - } else { - om <- matrix(0, R, T-1) } + else + om <- matrix(0, R, T-1) if(imm) { - iota <- predict(object, type="iota")[,1] + iota <- predict(object, type="iota", level=NULL, na.rm=FALSE)$Predicted iota <- matrix(iota, R, T-1, byrow=TRUE) - } else { - iota <- matrix(0, R, T-1) - } - srm <- object@sitesRemoved - if(length(srm) > 0){ - y <- y[-object@sitesRemoved,] - p <- p[-object@sitesRemoved,] # temporary workaround } + else + iota <- matrix(0, R, T-1) + ya <- array(y, c(R, J, T)) pa <- array(p, c(R, J, T)) post <- array(NA_real_, c(R, length(N), T)) @@ -947,7 +878,7 @@ setMethod("ranef", "unmarkedFitDailMadsen", } } for(i in 1:R) { - P <- matrix(1, K+1, K+1) + P <- matrix(1, K+1, K+1) switch(mix, P = g2 <- dpois(N, lam[i]), NB = { @@ -959,32 +890,14 @@ setMethod("ranef", "unmarkedFitDailMadsen", g2 <- (1-psi)*dpois(N, lam[i]) g2[1] <- psi + (1-psi)*exp(-lam[i]) }) - - #DETECTION MODEL - g1 <- rep(0, K+1) - cp <- pa[i,,1] - cp_na <- is.na(cp) - ysub <- ya[i,,1] - ysub[cp_na] <- NA - cp <- c(cp, 1-sum(cp, na.rm=TRUE)) - sumy <- sum(ysub, na.rm=TRUE) - - is_na <- c(is.na(ysub), FALSE) | is.na(cp) - - if(all(is.na(ysub))){ - post[i,,1] <- NA - } else { - - for(k in sumy:K){ - yit <- c(ysub, k-sumy) - g1[k+1] <- dmultinom(yit[!is_na], k, cp[!is_na]) - } - - g1g2 <- g1*g2 - post[i,,1] <- g1g2 / sum(g1g2) + g1 <- rep(1, K+1) + for(j in 1:J) { + if(is.na(ya[i,j,1]) | is.na(pa[i,j,1])) + next + g1 <- g1 * dbinom(ya[i,j,1], N, pa[i,j,1]) } - - + g1g2 <- g1*g2 + post[i,,1] <- g1g2 / sum(g1g2) for(t in 2:T) { if(!is.na(gam[i,t-1]) & !is.na(om[i,t-1])) { for(n0 in N) { @@ -1000,98 +913,20 @@ setMethod("ranef", "unmarkedFitDailMadsen", P <- P %*% P1 } } - - #DETECTION MODEL - g1 <- rep(0, K+1) - cp <- pa[i,,t] - cp_na <- is.na(cp) - ysub <- ya[i,,t] - ysub[cp_na] <- NA - cp <- c(cp, 1-sum(cp, na.rm=TRUE)) - sumy <- sum(ysub, na.rm=TRUE) - - is_na <- c(is.na(ysub), FALSE) | is.na(cp) - - if(all(is.na(ysub))){ - post[i,,t] <- NA - } else { - for(k in sumy:K){ - yit <- c(ysub, k-sumy) - g1[k+1] <- dmultinom(yit[!is_na], k, cp[!is_na]) - } - - g <- colSums(P * post[i,,t-1]) * g1 - post[i,,t] <- g / sum(g) + g1 <- rep(1, K+1) + for(j in 1:J) { + if(is.na(ya[i,j,t]) | is.na(pa[i,j,t])) + next + g1 <- g1 * dbinom(ya[i,j,t], N, pa[i,j,t]) } - } + g <- colSums(P * post[i,,t-1]) * g1 + post[i,,t] <- g / sum(g) + } } new("unmarkedRanef", post=post) }) -setMethod("ranef", "unmarkedFitNmixTTD", - function(object, ...) -{ - - M <- nrow(object@data@y) - J <- ncol(object@data@y) - tdist <- ifelse("shape" %in% names(object@estimates), "weibull", "exp") - mix <- ifelse("alpha" %in% names(object@estimates), "NB", "P") - K <- object@K - - yvec <- as.numeric(t(object@data@y)) - removed <- object@sitesRemoved - naflag <- as.numeric(is.na(yvec)) - surveyLength <- object@data@surveyLength - if(length(removed>0)) surveyLength <- surveyLength[-removed,] - ymax <- as.numeric(t(surveyLength)) - delta <- as.numeric(yvec