pcount: how to model siteCovs that vary across visits

[Corrado Alessandrini]

Dear all,

I want to model bird counts in line-transects, accounting for visit-specific detection variables (e.g. time of the survey, wind, cloud, disturbance..) and site-specific environmental variables (e.g. elevation, slope, solar radiation, NDVI..).

I carried out 3 visits at each transect (within late April and late June 2023), recording the abundance of my target species, let's say chaffinch.

The overall goal is to estimate the effect of my environemental covariates on the abundance of my target species, accounting for uneven detection probabilities, hence I think pcount() is what I need.

Now, the key point is that some siteCovs (namely, NDVI and solar radiation) DO VARY across the three visits: as the season progresses, they both increase at each site. To account for this, I calculated their values at the exact date of each visit (e.g. NDVI_1, NDVI_2, and NDVI_3 ..basically, as if they were obsCovs). 

And indeed they can easily be included as obsCovs, yet this would allow me to test their effect as detection variables (i.e. on the detection of chaffinches), rather than on their abundance (which is what I want).

So here is my question: How can I include NDVI and solar radiation as siteCovs in the unmarkedFrame?

Per each of these two variables, I guess that unmarkedFramePCount() should build a single column in the siteCovs section, to be used in the abundance formula in pcount(). (see script below)

Is there any way to do so?:)

Disclaimer: I'm new to unmarked (and this group too), so sorry if my question is silly!

Disclaimer2: I couldn't attach the script, so I paste it here below.

Any help or suggestion is very welcome!

Thanks a lot in advance and best regards,

Corrado Alessandrini

[ EXAMPLE SCRIPT ]

library(unmarked)
set.seed(16)

## 50 transects, visited 3 times each
df <- data.frame( transectID = c(1:50),
                 
                  n_birds_1 = rpois(50, lambda = 0.7), ## number of birds detected at each visit
                  n_birds_2 = rpois(50, lambda = 0.9),
                  n_birds_3 = rpois(50, lambda = 0.9),
                 
                  ## (visit-specific) detection variables    [obsCovs]
                  time_1  = rnorm(50, 6, 1),
                  time_2  = rnorm(50, 6, 1),
                  time_3  = rnorm(50, 6, 1),
                 
                  wind_1  = rpois(50, 1),
                  wind_2  = rpois(50, 1.2),
                  wind_3  = rpois(50, 2),
                 
                  ## (constant) site variables   [siteCovs]
                  elevation = rnorm(50, 1000, 300),
                  slope = rnorm(50, 45, 10),
                 
                  ## VISIT-SPECIFIC siteCovs <------------------ !!!
                  NDVI_1 = rnorm(50, 5000, 500),
                  NDVI_2 = rnorm(50, 6000, 700),
                  NDVI_3 = rnorm(50, 5500, 800)
                  )
                 
                 
df


## building unmarkedFrame

bird.udf <- with(df,{
  unmarkedFramePCount(
   
    y = cbind(n_birds_1,n_birds_2,n_birds_3), ## response variable
   
    ## (constant) site covs
    siteCovs = data.frame(elevation, slope
                          ,NDVI_1
                          ,NDVI_2
                          ,NDVI_3
                          ),
    ## (visit-specific) detection variables
    obsCovs = list(time = cbind(time_1,time_2,time_3),
                   wind = cbind(wind_1,wind_2,wind_3),
                   NDVI = cbind(NDVI_1,NDVI_2,NDVI_3)
                   )
  )})

summary(bird.udf)

## in siteCovs, NDVI is not a single variable,
## so, when modelling, it returns an error
mod1  = pcount( ~ scale(time) + wind
                ~ scale(elevation) + scale(slope) + scale(NDVI),
                mixture="P", K=30,
                bird.udf)
## "Error: Variable(s) NDVI not found in siteCovs"
#5: stop(paste("Variable(s)", paste(miss_vars, collapse = ", "),
#              "not found in siteCovs"), call. = FALSE)
#4: .local(umf, ...)
#3: getDesign(data, formula)
#2: getDesign(data, formula)
#1: pcount(~scale(time) + wind ~ scale(elevation) + scale(slope) +
#            scale(NDVI), mixture = "P", K = 30, bird.udf)


## adding NDVI as "cbind(NDVI_1,NDVI_2,NDVI_3)" in siteCovs still fails in producing a single column
bird.udf2 <- with(df,{
  unmarkedFramePCount(
   
    y = cbind(n_birds_1,n_birds_2,n_birds_3), ## response variable
   
    ## (constant) site covs
    siteCovs = data.frame(elevation, slope
                          ,NDVI = cbind(NDVI_1,NDVI_2,NDVI_3)
    ),
    ## (visit-specific) detection variables
    obsCovs = list(time = cbind(time_1,time_2,time_3),
                   wind = cbind(wind_1,wind_2,wind_3),
                   NDVI = cbind(NDVI_1,NDVI_2,NDVI_3)
    )
  )})
summary(bird.udf2) ## see, multiple columns for NDVI in siteCovs
## hence, same error is returned
mod1  = pcount( ~ scale(time) + wind
                ~ scale(elevation) + scale(slope) + scale(NDVI),
                mixture="P", K=30,
                bird.udf2)






## of course, placing NDVI in the detection formula DOES WORK,
##  but it is not what I need, since I want to model the effect of NDVI on the ABUNDANCE (not on the detection) of my target species
mod1  = pcount( ~ scale(time) + wind + scale(NDVI)
                ~ scale(elevation) + scale(slope),
                mixture="P", K=30,
                bird.udf)
summary(mod1)



## I think that using NDVI_1, NDVI_2, and NDVI_3 as three separate siteCovariates, e.g.
mod2  = pcount( ~ scale(time) + wind
                ~ scale(elevation) + scale(slope) + scale(NDVI_1) + scale(NDVI_2) + scale(NDVI_3),
                mixture="P", K=30,
                bird.udf)
summary(mod2)
## is not an option



## thanks a lot for going through this!!!

## CA - 11 June 2025

[Ken Kellner]

Hi Corrado,

This is not possible with a standard N-mixture model/pcount(). The model assumes constant abundance (and mean abundance, lambda) at a given site across all visits. If you had some kind of robust design where you had multiple primary periods between which NDVI varied, and then multiple repeated samples within each primary period, you could use the generalized N-mixture (gpcount) or open N-mixture (pcountOpen), I suppose, depending on your goals. But it doesn't sound like that would fit your design.

I'd recommend just averaging over the values of the covariates that vary to get a single value per site.

Ken

On Wed, Jun 11, 2025 at 07:41:54AM -0700, Corrado Alessandrini wrote:
> Dear all,
>
> I want to model bird counts in line-transects, accounting for
> visit-specific detection variables (e.g. time of the survey, wind, cloud,
> disturbance..) and site-specific environmental variables (e.g. elevation,
> slope, solar radiation, NDVI..).
> I carried out 3 visits at each transect (within late April and late June
> 2023), recording the abundance of my target species, let's say chaffinch.
>
> The overall goal is to estimate the effect of my environemental covariates
> on the abundance of my target species, accounting for uneven detection
> probabilities, hence I think pcount() is what I need.
>

> Now, the key point is that *some siteCovs (namely, NDVI and solar
> radiation) DO VARY across the three visits*: as the season progresses, they

> both increase at each site. To account for this, I calculated their values
> at the exact date of each visit (e.g. NDVI_1, NDVI_2, and NDVI_3
> ..basically, as if they were obsCovs).
> And indeed they can easily be included as obsCovs, yet this would allow me
> to test their effect as detection variables (i.e. on the detection of
> chaffinches), rather than on their abundance (which is what I want).
> So here is my question: How can I include NDVI and solar radiation as
> siteCovs in the unmarkedFrame?

> Per each of these two variables, I guess that *unmarkedFramePCount()*

> should build a single column in the siteCovs section, to be used in the
> abundance formula in pcount(). (see script below)
> Is there any way to do so?:)
>
> Disclaimer: I'm new to unmarked (and this group too), so sorry if my
> question is silly!
> Disclaimer2: I couldn't attach the script, so I paste it here below.
>
> Any help or suggestion is very welcome!
> Thanks a lot in advance and best regards,
> Corrado Alessandrini
>
>

> *[ EXAMPLE SCRIPT ]*

> --
> *** Three hierarchical modeling email lists ***
> (1) unmarked (this list): for questions specific to the R package unmarked
> (2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
> (3) HMecology: for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
> ---
> You received this message because you are subscribed to the Google Groups "unmarked" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/unmarked/b0f618a3-0be9-422c-8d69-5070fa25c3a1n%40googlegroups.com.

[Marc Kery]

Dear Corrado,

the short answer is this: get the average of these covariates over the 3 visits and use the resulting vectors as the site level covariates.

The longer answer is this: inherent in the use of the N-mixture model for a single season is the closure assumption: i.e., that abundance at site i, N_i, does not change over repeated visits. Therefore, it would not make sense to try and use a time-varying covariate to explain variation in N, since we assume that N can only vary in space, but not in time. Therefore, for the simplest, single-season Nmix model, the covariates for N can be only be spatial, while covariates for p may vary both in space (among sites) and in time (among visits).

Best regards  ----- Marc

---

From: unma...@googlegroups.com <unma...@googlegroups.com> on behalf of Corrado Alessandrini <corrado.al...@herbando.com>
Sent: Wednesday, June 11, 2025 16:41
To: unmarked <unma...@googlegroups.com>
Subject: [unmarked] pcount: how to model siteCovs that vary across visits

 

[Corrado Alessandrini]

Dear Ken and Marc,

thanks a lot for your cristal-clear answer.

I now defenitely see the point! The closure assumption would be necessarily violated if I look for a time-varing effect of any of my siteCovs. The sampling design proposed by Ken would overcome this issue, but sadly I'm not in that situation. I will eventually go for the averaging option :)

Thanks a lot again!

Bests,

Corrado

[Anam Ahsan]

Hi 

Please can anyone help me 

Like I have a top model with lowest aic  for tiger occupancy 

Now I find that covariates forest and water is the main predictors in occupancy 

Now I want to find one psi for tiger 

As well as predicted psi value at site covariates forest and water 

How can I find the one psi for tiger like is it the mean of all psi in all my grids ?

Anam 🍀

On Jun 12, 2025, at 05:16, Corrado Alessandrini <corrado.al...@herbando.com> wrote:

Dear Ken and Marc,

To view this discussion visit https://groups.google.com/d/msgid/unmarked/510a7d65-0a1b-4f85-8012-c9da493af289n%40googlegroups.com.