Skip to content

Commit

Permalink
Add tests and fix param ordering bug
Browse files Browse the repository at this point in the history
  • Loading branch information
@kenkellner
kenkellner committed Jan 28, 2025
1 parent db2853d commit 13eca4a
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 27 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -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(
Expand Down
50 changes: 25 additions & 25 deletions R/IDS.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)){
Expand Down Expand Up @@ -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])
Expand All @@ -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)){
Expand All @@ -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){
Expand All @@ -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(),
Expand Down
86 changes: 86 additions & 0 deletions tests/testthat/test_IDS.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
})

0 comments on commit 13eca4a

Please sign in to comment.