Practical 3: Dealing with false positives
Olivier Gimenez (adapted from David Miller’s code)
May 2023
1 Objective
- Fit (static) occupancy models accomodating false positives
2 Pre-requisites
Load some packages:
library(unmarked)3 Data
3.1 Description
We use data from the paper Determining Occurrence Dynamics when False Positives Occur: Estimating the Range Dynamics of Wolves from Public Survey Data by Miller et al. published in Plos One and available at https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0065808.
Data are of 2 types. There are data from hunters who report whether they observed wolves while hunting for deer or elk. Other data come from observations made by state agency personnel based on collaring and relocation by radio telemetry. Covariates include the amount of hunter effort, habitat quality, and whether observations were by hunters or agency personnel. Static analysis uses data from the year 2010.
We thank David Miller for sharing the data.
In the figure below, we display the data. Cells were divided based on
perceived habitat quality into low, medium, and high quality categories
(A). Most certain observations of known packs made by state agency
personnel and based on collaring and relocation by radio telemetry were
concentrated in the western end of the study area where the higher
quality habitat was located (B). While high frequency of observations by
hunters also occurred in high-quality areas, they also reported a low
frequency of wolf observations in the eastern portion of the study area,
which we suspected were due to mis-identification (C).
3.2 Read in data
Read in data:
wolves <- read.csv("data/MTwolves.csv")Inspect data:
str(wolves)## 'data.frame': 257 obs. of 52 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ X2007NA1: int 0 0 1 0 1 0 0 0 0 0 ...
## $ X2007NA2: int 0 0 0 0 0 1 0 0 0 0 ...
## $ X2007NA3: int 0 1 0 0 1 1 0 0 0 1 ...
## $ X2007NA4: int 0 0 1 0 0 1 1 0 0 1 ...
## $ X2007NA5: int 0 0 0 1 1 0 1 0 0 1 ...
## $ X2007NAK: int 0 0 0 0 0 0 0 0 0 2 ...
## $ X2008NA1: int 0 0 1 1 1 1 1 0 0 1 ...
## $ X2008NA2: int 0 0 1 0 1 1 0 1 0 1 ...
## $ X2008NA3: int 0 0 0 1 1 1 1 0 0 1 ...
## $ X2008NA4: int 0 0 0 1 1 1 1 0 0 1 ...
## $ X2008NA5: int 0 1 1 1 1 1 1 0 1 1 ...
## $ X2008NAK: int 0 0 0 0 2 0 2 0 0 0 ...
## $ X2009NA1: int 0 0 1 1 1 1 0 0 0 1 ...
## $ X2009NA2: int 0 0 1 1 1 1 1 0 0 1 ...
## $ X2009NA3: int 1 1 1 1 1 1 1 0 0 1 ...
## $ X2009NA4: int 0 0 1 1 1 0 0 0 0 1 ...
## $ X2009NA5: int 0 0 0 0 1 1 1 0 0 1 ...
## $ X2009NAK: int 0 0 0 0 2 0 2 0 2 2 ...
## $ X2010NA1: int 0 0 0 1 1 1 1 0 0 1 ...
## $ X2010NA2: int 0 0 1 1 1 1 1 1 0 1 ...
## $ X2010NA3: int 0 1 1 1 1 1 1 1 0 1 ...
## $ X2010NA4: int 0 0 0 1 1 1 1 0 0 1 ...
## $ X2010NA5: int 0 0 1 1 1 1 1 0 0 1 ...
## $ X2010NAK: int 0 0 0 0 2 2 2 0 2 2 ...
## $ High : int 1 1 1 1 1 0 1 1 1 1 ...
## $ Medium : int 0 0 0 0 0 1 0 0 0 0 ...
## $ Low : int 0 0 0 0 0 0 0 0 0 0 ...
## $ E07NA1 : num 14.4 21.9 22 21.4 16.8 ...
## $ E07NA2 : num 14.4 21.9 22 21.4 16.8 ...
## $ E07NA3 : num 14.4 21.9 22 21.4 16.8 ...
## $ E07NA4 : num 14.4 21.9 22 21.4 16.8 ...
## $ E07NA5 : num 14.4 21.9 22 21.4 16.8 ...
## $ E07NAK : num 0 0 0 0 0 0 0 0 0 0 ...
## $ E08NA1 : num 14.9 22 22.1 21.5 16.8 ...
## $ E08NA2 : num 14.9 22 22.1 21.5 16.8 ...
## $ E08NA3 : num 14.9 22 22.1 21.5 16.8 ...
## $ E08NA4 : num 14.9 22 22.1 21.5 16.8 ...
## $ E08NA5 : num 14.9 22 22.1 21.5 16.8 ...
## $ E08NAK : num 0 0 0 0 0 0 0 0 0 0 ...
## $ E09NA1 : num 16.2 20.5 20.5 19.7 13.7 ...
## $ E09NA2 : num 16.2 20.5 20.5 19.7 13.7 ...
## $ E09NA3 : num 16.2 20.5 20.5 19.7 13.7 ...
## $ E09NA4 : num 16.2 20.5 20.5 19.7 13.7 ...
## $ E09NA5 : num 16.2 20.5 20.5 19.7 13.7 ...
## $ E09NAK : num 0 0 0 0 0 0 0 0 0 0 ...
## $ E10NA1 : num 13.5 20.8 20.9 19.9 12.7 ...
## $ E10NA2 : num 13.5 20.8 20.9 19.9 12.7 ...
## $ E10NA3 : num 13.5 20.8 20.9 19.9 12.7 ...
## $ E10NA4 : num 13.5 20.8 20.9 19.9 12.7 ...
## $ E10NA5 : num 13.5 20.8 20.9 19.9 12.7 ...
## $ E10NAK : num 0 0 0 0 0 0 0 0 0 0 ...4 Modelling
4.1 Start with only hunter reporting data for wolves and no accounting for false positives
Build list of covariate:
oc <- list(effort = wolves[,47:51])Organise data and covariates:
umf <- unmarkedFrameOccu(y = wolves[,20:24], obsCovs = oc)Summarize data:
str(umf)## Formal class 'unmarkedFrameOccu' [package "unmarked"] with 5 slots
## ..@ y : int [1:257, 1:5] 0 0 0 1 1 1 1 0 0 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:5] "X2010NA1" "X2010NA2" "X2010NA3" "X2010NA4" ...
## ..@ obsCovs :'data.frame': 1285 obs. of 1 variable:
## .. ..$ effort: num [1:1285] 13.5 13.5 13.5 13.5 13.5 ...
## ..@ siteCovs: NULL
## ..@ mapInfo : NULL
## ..@ obsToY : num [1:5, 1:5] 1 0 0 0 0 0 1 0 0 0 ...Fit model:
fm1 <- occu(~effort ~ 1, data = umf)Display model coefficients (parameters)
coef(fm1)## psi(Int) p(Int) p(effort)
## 2.9393910 -2.5989479 0.3643864Get occupancy probability estimate:
1/(1 + exp(-coef(fm1)[1]))## psi(Int)
## 0.9497597Or:
plogis(coef(fm1)[1])## psi(Int)
## 0.9497597Or:
backTransform(fm1, 'state')## Backtransformed linear combination(s) of Occupancy estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.95 0.0354 2.94 1
##
## Transformation: logisticNow let’s plot the detection probability for hunters as a function of effort:
effort_grid <- seq(min(oc$effort), max(oc$effort), length = 100) # grid
logitp <- coef(fm1)[2] + coef(fm1)[3] * effort_grid # predictions on logit scale
p <- plogis(logitp) # back-transformation
plot(effort_grid,
p,
type = "l",
lwd = 3,
xlab = "effort",
ylab = "detection probability",
ylim = c(0,1))
Compare figure below with Figure 3 in Miller and colleagues’ paper.
4.2 Still only data for hunters but account for false positives
Get covariates:
effort <- matrix(unlist(wolves[,47:51]), nrow(wolves), 5)
oc <- list(effort = effort)Organise data and covariates using unmarkedFrameOccuFP
function. Note the extra argument type which takes a vector
with 3 values designating the number of occasions where data is of type
1, type 2, and type 3. For type 1 data, false positive detections are
assumed not to occur. For type 2 data, both false negative and false
positive detection probabilities are estimated. Both p or p11 (the true
positive detection probability) and fp or p01 (the false positive
detection probability) are estimated for occasions when this data type
occurs. For type 3 data, observations are assumed to include both
certain detections (false positives assumed not to occur) and uncertain
detections that may include false positive detections. Let’s get to
it:
umfFP <- unmarkedFrameOccuFP(y = wolves[,20:24],
obsCovs = oc,
type = c(0,5,0)) # 5 columns are of type 2Summarize data:
str(umfFP) ## Formal class 'unmarkedFrameOccuFP' [package "unmarked"] with 6 slots
## ..@ type : num [1:3] 0 5 0
## ..@ y : int [1:257, 1:5] 0 0 0 1 1 1 1 0 0 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:5] "X2010NA1" "X2010NA2" "X2010NA3" "X2010NA4" ...
## ..@ obsCovs :'data.frame': 1285 obs. of 1 variable:
## .. ..$ effort: num [1:1285] 13.5 13.5 13.5 13.5 13.5 ...
## ..@ siteCovs: NULL
## ..@ mapInfo : NULL
## ..@ obsToY : num [1:5, 1:5] 1 0 0 0 0 0 1 0 0 0 ...Fit model with effort on both detection and false positive detection probabilities:
fm2 <- occuFP(detformula = ~effort, # detection
FPformula = ~effort, # false positive detection
Bformula = ~1, # probability detections are certain
stateformula = ~1, # occupancy
starts = c(0, 0, 0, -5, 0), # initial values
data = umfFP)Inspect estimates:
coef(fm2)## psi(Int) p(Int) p(effort) fp(Int) fp(effort)
## -0.7651619 -1.4539025 0.1380600 -7.4198442 1.7210737Get occupancy prob estimate:
backTransform(fm2, 'state')## Backtransformed linear combination(s) of Occupancy estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.318 0.0621 -0.765 1
##
## Transformation: logistic4.3 Add in the agency data
Build covariates:
effort <- matrix(unlist(wolves[,c(52,47:51)]), nrow(wolves), 6)
oc <- list(effort = effort)Get detections and non-detections:
y <- wolves[,c(25,20:24)]
y[,1] <- y[,1] / 2Organise data and covariates:
umfFP <- unmarkedFrameOccuFP(y = y,
obsCovs = oc,
type = c(1,5,0)) # first column of type 1, 5 columns of type 2Summarize data:
str(umfFP) ## Formal class 'unmarkedFrameOccuFP' [package "unmarked"] with 6 slots
## ..@ type : num [1:3] 1 5 0
## ..@ y : num [1:257, 1:6] 0 0 0 0 1 1 1 0 1 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:6] "X2010NAK" "X2010NA1" "X2010NA2" "X2010NA3" ...
## ..@ obsCovs :'data.frame': 1542 obs. of 1 variable:
## .. ..$ effort: num [1:1542] 0 13.5 13.5 13.5 13.5 ...
## ..@ siteCovs: NULL
## ..@ mapInfo : NULL
## ..@ obsToY : num [1:6, 1:6] 1 0 0 0 0 0 0 1 0 0 ...Fit model:
fm3 <- occuFP(detformula = ~effort, # detection
FPformula = ~effort, # false positives detection
Bformula = ~1, # probability detections are certain
stateformula = ~1, # occupancy
starts = c(0, 0, 0, -5, 0), # initial values
data = umfFP) Print results:
fm3##
## Call:
## occuFP(detformula = ~effort, FPformula = ~effort, Bformula = ~1,
## stateformula = ~1, data = umfFP, starts = c(0, 0, 0, -5,
## 0))
##
## Occupancy:
## Estimate SE z P(>|z|)
## -0.63 0.187 -3.37 0.000742
##
## Detection:
## Estimate SE z P(>|z|)
## (Intercept) -0.126 0.2700 -0.468 6.40e-01
## effort 0.205 0.0304 6.761 1.37e-11
##
## false positive:
## Estimate SE z P(>|z|)
## (Intercept) -3.520 0.2385 -14.8 2.58e-49
## effort 0.303 0.0287 10.6 3.54e-26
##
## AIC: 1146.8764.4 Account for habitat quality
Get covariates:
effort <- matrix(unlist(wolves[,c(52,47:51)]), nrow(wolves), 6)
oc <- list(effort = effort)Get detections and non-detections:
y <- wolves[,c(25,20:24)]
y[,1] <- y[,1] / 2Organise data and covariates:
umfFP <- unmarkedFrameOccuFP(y = y,
siteCovs = wolves[,26:28],
obsCovs = oc,
type = c(1,5,0)) # first column of type 1, 5 columns of type 2Summarize data:
str(umfFP) ## Formal class 'unmarkedFrameOccuFP' [package "unmarked"] with 6 slots
## ..@ type : num [1:3] 1 5 0
## ..@ y : num [1:257, 1:6] 0 0 0 0 1 1 1 0 1 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:6] "X2010NAK" "X2010NA1" "X2010NA2" "X2010NA3" ...
## ..@ obsCovs :'data.frame': 1542 obs. of 1 variable:
## .. ..$ effort: num [1:1542] 0 13.5 13.5 13.5 13.5 ...
## ..@ siteCovs:'data.frame': 257 obs. of 3 variables:
## .. ..$ High : int [1:257] 1 1 1 1 1 0 1 1 1 1 ...
## .. ..$ Medium: int [1:257] 0 0 0 0 0 1 0 0 0 0 ...
## .. ..$ Low : int [1:257] 0 0 0 0 0 0 0 0 0 0 ...
## ..@ mapInfo : NULL
## ..@ obsToY : num [1:6, 1:6] 1 0 0 0 0 0 0 1 0 0 ...Fit model:
fm4 <- occuFP(detformula = ~effort, # detection
FPformula = ~effort, # false positives detection
Bformula = ~1, # probability detections are certain
stateformula = ~High + Medium + Low - 1, # occupancy
starts = c(0, 0, 0, 0, 0, -5, 0), # initial values
data = umfFP)Occupancy estimates for each level of habitat quality, starting by high:
lc <- linearComb(fm4, c(1, 0, 0), type="state") # occ on the logit scale when high
psi_high <- backTransform(lc)Then medium:
lc <- linearComb(fm4, c(0, 1, 0), type="state") # occ on the logit scale when medium
psi_medium <- backTransform(lc)Eventually low:
lc <- linearComb(fm4, c(0, 0, 1), type="state") # occ on the logit scale when low
psi_low <- backTransform(lc)Compare the following estimates with Figure 4, upper left panel:
psi_high@estimate## [1] 0.8347993psi_medium@estimate## [1] 0.9040658psi_low@estimate## [1] 0.03511551Get detection estimates as a function of effort:
effort_grid <- seq(min(oc$effort), max(oc$effort), length = 100)
p11 <- plogis(coef(fm4)[4] + coef(fm4)[5] * effort_grid) # detection
p01 <- plogis(coef(fm4)[6] + coef(fm4)[7] * effort_grid) # false positives detectionCompare figure below with Figure 3:
plot(effort_grid,
p11,
type = "l",
lwd = 3,
xlab = "effort",
ylab = "detection probability",
col = "red",
ylim = c(0,1))
lines(effort_grid,
p01,
lwd = 3,
col = "blue")
legend("bottomright", legend = c("true positive det prob (p11)",
"false positive det prob (p01)"),
col = c("red", "blue"),
lwd = 3)
Both the true positive (p11) and false positive (p01) detection probabilities increase as hunter effort increases (Figure 3). False positive probabilities are greater than 0.5 in areas with the greatest hunting effort.