diff --git a/DESCRIPTION b/DESCRIPTION index d0fe0c2e..26140b1a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.4.3.9007 -Date: 2025-01-09 +Version: 1.4.3.9008 +Date: 2025-01-28 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( diff --git a/R/IDS.R b/R/IDS.R index 6cb6edc1..e97510d9 100644 --- a/R/IDS.R +++ b/R/IDS.R @@ -168,7 +168,7 @@ IDS <- function(lambdaformula = ~1, # Hazard-rate scale parameter(s) if(keyfun == "hazard"){ pind_mat[6,] <- max(pind_mat) + 1 - if(!is.null(detformulaOC) & !is.null(dataOC)){ + if(!is.null(detformulaPC) & !is.null(dataPC)){ pind_mat[7,] <- max(pind_mat) + 1 } if(!is.null(detformulaOC) & !is.null(dataOC)){ @@ -278,14 +278,6 @@ IDS <- function(lambdaformula = ~1, est_list <- list(lam=lam_est, ds=dist_est) - if(keyfun == "hazard"){ - ds_haz_coef <- get_coef_info(sdr, "schds", "(Intercept)", pind_mat[6,1]:pind_mat[6,2]) - ds_est_haz <- unmarkedEstimate("Distance sampling scale", short.name="ds_scale", - estimates = ds_haz_coef$ests[1], covMat = ds_haz_coef$cov, fixed=1, - invlink="exp", invlinkGrad="exp") - est_list <- c(est_list, list(schds=ds_est_haz)) - } - if(!is.null(detformulaPC) & !is.null(dataPC)){ pc_coef <- get_coef_info(sdr, "pc", colnames(gd_pc$V), pind_mat[3,1]:pind_mat[3,2]) @@ -294,14 +286,6 @@ IDS <- function(lambdaformula = ~1, estimates = pc_coef$ests, covMat = pc_coef$cov, fixed=1:ncol(gd_pc$V), invlink = "exp", invlinkGrad = "exp") est_list <- c(est_list, list(pc=pc_est)) - - if(keyfun == "hazard"){ - pc_haz_coef <- get_coef_info(sdr, "scpc", "(Intercept)", pind_mat[7,1]:pind_mat[7,2]) - pc_est_haz <- unmarkedEstimate("Point count scale", short.name="pc_scale", - estimates = pc_haz_coef$ests[1], covMat = pc_haz_coef$cov, fixed=1, - invlink="exp", invlinkGrad="exp") - est_list <- c(est_list, list(scpc=pc_est_haz)) - } } if(!is.null(detformulaOC) & !is.null(dataOC)){ @@ -311,14 +295,6 @@ IDS <- function(lambdaformula = ~1, estimates = oc_coef$ests, covMat = oc_coef$cov, fixed=1:ncol(gd_oc$V), invlink = "exp", invlinkGrad = "exp") est_list <- c(est_list, list(oc=oc_est)) - - if(keyfun == "hazard"){ - oc_haz_coef <- get_coef_info(sdr, "scoc", "(Intercept)", pind_mat[8,1]:pind_mat[8,2]) - oc_est_haz <- unmarkedEstimate("Presence/absence scale", short.name="oc_scale", - estimates = oc_haz_coef$ests[1], covMat = oc_haz_coef$cov, fixed=1, - invlink="exp", invlinkGrad="exp") - est_list <- c(est_list, list(scoc=oc_est_haz)) - } } if(has_avail){ @@ -330,6 +306,30 @@ IDS <- function(lambdaformula = ~1, est_list <- c(est_list, list(phi=avail_est)) } + if(keyfun == "hazard"){ + ds_haz_coef <- get_coef_info(sdr, "schds", "(Intercept)", pind_mat[6,1]:pind_mat[6,2]) + ds_est_haz <- unmarkedEstimate("Distance sampling scale", short.name="ds_scale", + estimates = ds_haz_coef$ests[1], covMat = ds_haz_coef$cov, fixed=1, + invlink="exp", invlinkGrad="exp") + est_list <- c(est_list, list(schds=ds_est_haz)) + + if(!is.null(detformulaPC) & !is.null(dataPC)){ + pc_haz_coef <- get_coef_info(sdr, "scpc", "(Intercept)", pind_mat[7,1]:pind_mat[7,2]) + pc_est_haz <- unmarkedEstimate("Point count scale", short.name="pc_scale", + estimates = pc_haz_coef$ests[1], covMat = pc_haz_coef$cov, fixed=1, + invlink="exp", invlinkGrad="exp") + est_list <- c(est_list, list(scpc=pc_est_haz)) + } + + if(!is.null(detformulaOC) & !is.null(dataOC)){ + oc_haz_coef <- get_coef_info(sdr, "scoc", "(Intercept)", pind_mat[8,1]:pind_mat[8,2]) + oc_est_haz <- unmarkedEstimate("Presence/absence scale", short.name="oc_scale", + estimates = oc_haz_coef$ests[1], covMat = oc_haz_coef$cov, fixed=1, + invlink="exp", invlinkGrad="exp") + est_list <- c(est_list, list(scoc=oc_est_haz)) + } + } + est_list <- unmarkedEstimateList(est_list) new("unmarkedFitIDS", fitType = "IDS", call = match.call(), diff --git a/tests/testthat/test_IDS.R b/tests/testthat/test_IDS.R index f9281ee8..fa7b0a64 100644 --- a/tests/testthat/test_IDS.R +++ b/tests/testthat/test_IDS.R @@ -206,3 +206,89 @@ test_that("IDS handles missing values", { )) }) + +test_that("Hazard-rate works with IDS", { + + # Simulate IDS dataset + set.seed(123) + + # First simulate distance sampling data + # Distance sampling sites + Mds <- 300 + elev <- rnorm(Mds) + # Distance breaks + db <- seq(0, 0.3, length.out=7) + # blank unmarked frame + umf_ds <- unmarkedFrameDS(y = matrix(NA, Mds, length(db)-1), + siteCovs = data.frame(elev = elev), + survey="point", dist.breaks = db, unitsIn='km') + + # True coefficient values + cf <- list(state = c(3, -0.5), # abundance intercept and effect of elevation + det = -2.5, # detection sigma - exp(-2.5) = 0.08 + scale = log(2)) # detection hazard rate scale = 2 + + # Simulate response + umf_ds <- simulate(umf_ds, formula=~1~elev, coefs=cf, + keyfun = "hazard", output = "density", unitsOut = "kmsq")[[1]] + + # Fit regular distance sampling model as test + mod <- distsamp(~1~elev, umf_ds, keyfun="hazard", output="density", unitsOut = "kmsq") + mod + + # good correspondance + expect_equivalent(coef(mod), c(2.8913, -0.4381, -2.4697, 0.60489), tol=1e-4) + + # "Point count" sites + # simulate these as 1-bin distance + Mpc <- 500 + elev2 <- rnorm(Mpc) + max_dist <- 0.5 + db2 <- c(0, max_dist) # single bin + + # blank unmarked frame + umf_pc <- unmarkedFrameDS(y = matrix(NA, Mpc, 1), # single column in y + siteCovs = data.frame(elev = elev2), + survey = "point", dist.breaks=db2, unitsIn="km") + + # simulate response (same coefficients) + umf_pc <- simulate(umf_pc, formula=~1~elev, coefs=cf, + keyfun = "hazard", output = "density", unitsOut = "kmsq")[[1]] + + + # convert to unmarkedFramePcount + umf_pc2 <- unmarkedFramePCount(y = umf_pc@y, + siteCovs = umf_pc@siteCovs) + + # Fit IDS model (in this case: common abundance and detection processes, + # and common duration) + mod2 <- IDS(lambdaformula = ~elev, detformulaDS = ~1, + dataDS = umf_ds, dataPC = umf_pc2, keyfun = "hazard", + unitsOut = "kmsq", maxDistPC = 0.5) + + # similar to just distance sampling data + expect_equivalent(coef(mod2), + c(2.8419,-0.5070, -2.4420, 0.6313), tol=1e-5) + #cbind(true=unlist(cf), est=coef(mod2)) + + # Check with separate formulas + mod3 <- IDS(lambdaformula = ~elev, detformulaDS = ~1, detformulaPC = ~1, + dataDS = umf_ds, dataPC = umf_pc2, keyfun = "hazard", + unitsOut = "kmsq", maxDistPC = 0.5) + + # Also make sure the coefs are in the right order: lam, then all sigma, then all haz scales + expect_equal(coef(mod3), + c(`lam(Int)` = 2.86807, `lam(elev)` = -0.50746, + `ds(Int)` = -2.47106, `pc(Int)` = -3.37866, + `ds_scale(Int)` = 0.60447, `pc_scale(Int)` = -0.02817), tol=1e-5) + + # estimation doesn't really work for this one + expect_warning(se <- SE(mod3)) + expect_true(is.nan(se["pc(Int)"])) + expect_true(is.nan(se["pc_scale(Int)"])) + + # Make sure different formulas for sigma work + mod4 <- IDS(lambdaformula = ~1, detformulaDS = ~elev, detformulaPC = ~1, + dataDS = umf_ds, dataPC = umf_pc2, keyfun = "hazard", + unitsOut = "kmsq", maxDistPC = 0.5) +})