Trend models in pcountOpen
fle...@ag.arizona.edu
Dear Unmarked Community:
I am hoping to use pcountOpen to estimate 1) temporal trends in abundance and 2) the effects annual weather factors on temporal variation in abundance of various lizard species in the Sonoran Desert. I am not great with R and fairly new to Unmarked so my questions should be representative of many prospective users (and hopefully not too naïve!). After reading all posts with the search term “pcountOpen” and the webinar/workshop materials, I have fit trend models that seem to be providing reasonable estimates but have questions on model fitting, strategy, and interpretation that I would be very grateful for feedback on.
Background: Lizards were surveyed along up to 32 fixed transects per year during one daily visit a year in both spring and summer for 25 years. Sixteen transects were surveyed most years, and 16 others were added later in the study (mean = 23 transects/yr) and sampled annually depending on resources. During each daily visit, transects were surveyed repeatedly 3-8 times from ~6AM to 12PM to ensure overlap with times of peak animal activity. Thus, I assumed closure among repeated surveys during each daily visit. Because closure cannot be reliably assumed between seasons in the same year, I developed a separate model for each species in each season. Summary of an unmarked frame for one species and season, and the best detection model follows:
32 sitesMaximum number of observations per site: 200
Mean number of observations per site: 84.47
Number of primary survey periods: 25
Number of secondary survey periods: 8
Sites with at least one detection: 32
Tabulation of y observations:
0 1 2 3 4 5 6 7 8 9 10 11 12 14 15 <NA>
1234 647 373 198 111 61 40 12 15 4 3 1 2 1 1 3697
Site-level covariates:
Hydro Topo Soil Subdiv Veg
mixed : 3 bajada : 8 coarse:12 AZUP : 8 shrub :10
upland :16 rocky slope : 7 fine :20 LCRV :14 woodland:22
xeroriparian:13 valley floor:17 Trans:10
Observation-level covariates:
Temp MRain DOY TAS
Min. :-2.358 Min. :-0.951 Min. :-2.642 Min. :-2.026
1st Qu.:-0.805 1st Qu.:-0.834 1st Qu.:-0.800 1st Qu.:-0.830
Median :-0.038 Median :-0.287 Median : 0.149 Median :-0.006
Mean : 0.000 Mean : 0.000 Mean : 0.000 Mean : 0.000
3rd Qu.: 0.823 3rd Qu.: 0.613 3rd Qu.: 0.819 3rd Qu.: 0.781
Max. : 2.273 Max. : 4.859 Max. : 2.158 Max. : 3.786
NA's :3697 NA's :3697 NA's :3697 NA's :3697
Yearly-site-level covariates:
Yr YearNom
Min. :-12 Yr00 : 32
1st Qu.: -6 Yr01 : 32
Median : 0 Yr02 : 32
Mean : 0 Yr03 : 32
3rd Qu.: 6 Yr04 : 32
Max. : 12 Yr05 : 32
(Other):608
BEST DETECTION MODEL
Call:
pcountOpen(lambdaformula = ~1, gammaformula = ~1, omegaformula = ~1, pformula = ~Temp + I(Temp^2) + TAS + I(TAS^2) + Topo + MRain,
data = cnti.umf, mixture = "P", K = 60, dynamics = "trend", se = T)
Abundance:
Estimate SE z P(>|z|)
2.61 0.126 20.7 3.2e-95
Recruitment:
Estimate SE z P(>|z|)
-0.00484 0.0115 -0.42 0.675
Detection:
Estimate SE z P(>|z|)
(Intercept) -2.1559 0.1698 -12.70 6.26e-37
Temp -0.0979 0.0575 -1.70 8.84e-02
I(Temp^2) -0.3140 0.0307 -10.24 1.26e-24
TAS 0.3547 0.0562 6.31 2.78e-10
I(TAS^2) -0.0657 0.0262 -2.51 1.21e-02
Toporocky slope -0.8990 0.2029 -4.43 9.43e-06
Topovalley floor 0.2677 0.1550 1.73 8.42e-02
MRain -0.0722 0.0245 -2.94 3.24e-03
1) Is there a better way to handle data from each season? Because the presence and magnitude of trends vary between seasons for the same species, I’d like assess trends for both seasons together but am unsure how to proceed in one model. Including a row for each transect in each season violates the independence assumption, while pooling data among secondary periods for both seasons in each year violates the closure assumption. Is there a way to do this or is it best to model each season separately?
2) Is it necessary to model the initial abundance (lambda) when fitting trend models (dynamics=“trend”). One somewhat related post suggests the need to fit a siteCov with values from year 1 in lambda formula, while fitting a yearlySiteCovs for all years in gamma formula (e.g., pcountOpen(lambdaformula = ~Yr1, gammaformula = ~Yr1_25, …). In my case, I input yearlySiteCovs for year as 25 columns for each year with dummy variables ranging from -12 to 12 centered on the middle year=0.
3) The manual indicates the trend parameter gamma is the finite rate of population change (e.g., lambda), but my output suggests it is the instantaneous rate of population change (r). For example, the trend model below estimates a significant gamma of -0.0123. Thus, based on the calculations below, I assume lambda = 0.9878 and that populations are declining at a rate = 1.22%/yr or 25.5% over all years. Which interpretation is correct?
> (fm_best <- pcountOpen(~1, ~Yr, ~1, ~Temp + I(Temp^2) + TAS + I(TAS^2)+ Topo + MRain , cnti.umf, K=70, mixture = "P", dynamics = "trend", se = TRUE))
Call:
pcountOpen(lambdaformula = ~1, gammaformula = ~Yr, omegaformula = ~1, pformula = ~Temp + I(Temp^2) + TAS + I(TAS^2) + Topo + MRain,
data = cnti.umf, mixture = "P", K = 70, dynamics = "trend", se = TRUE)
Abundance:
Estimate SE z P(>|z|)
2.47 0.21 11.8 4.17e-32
Recruitment:
Estimate SE z P(>|z|)
(Intercept) -0.0505 0.02457 -2.06 0.03975
Yr -0.0123 0.00468 -2.63 0.00856
Detection:
Estimate SE z P(>|z|)
(Intercept) -2.9887 0.4656 -6.419 1.37e-10
Temp 0.0995 0.1054 0.944 3.45e-01
I(Temp^2) -0.3329 0.0528 -6.300 2.97e-10
TAS 0.2351 0.1019 2.306 2.11e-02
I(TAS^2) -0.0787 0.0462 -1.704 8.84e-02
Toporocky slope 0.5512 0.4839 1.139 2.55e-01
Topovalley floor 1.0928 0.4503 2.427 1.52e-02
MRain -0.0733 0.0405 -1.811 7.02e-02
AIC: 28613.12
LambdaSum <- (exp(-0.0123)) # lambda = 0.9878 or a 1.22% decline/yr
Total_declineSu <- (1-(0.9878^24))*100 # 25.52% decline
4) My final question concerns modifying my trend models to address the 2nd objective (how weather affects temporal variation in abundance). To do so, do I maintain the dynamics=“trend” syntax and simply replace my yearlySiteCovs (and perhaps SiteCovs for initial abundance) for year with those for weather factors? I have explored other fitting options (e.g., dynamics=“autoreg”) but they run for many hrs and then typically crash R. If anyone can provide some example code or suggestions for doing this I would be very grateful.
Thanks tons!
Aaron
Jeff
I'm not sure what the best approach is for the two seasons, but I may be able to help with a few of the other questions.
2) I don't think this is necessary. Any covariates for initial abundance/lambda should be ones that differ between sites in the initial year and could affect abundance. Fitting a year effect to initial abundance doesn't really make sense, because all sites in year 1 will have the same value(s) for this.
3) I don't think either interpretation is quite correct. The trend is indeed on the log scale, and should be exponentiated to get finite rate of change (in most places known as lambda). But the model you show for this question does not give the overall population growth rate on either scale; rather, you're modeling a population growth rate that is itself declining by year (and that this trend is significant). The first model you showed gives the overall population growth rate (r = -0.00484, lambda = 0.99517), which is not significantly different from stability. From your description, it looks like the second model has the following population growth rates:
year r
1 -12 0.0971 (-0.0505 + -12*-0.0123)
2 -11 0.0848 (-0.0505 + -11*-0.0123)
3 -10 0.0725
4 -9 0.0602
5 -8 0.0479
6 -7 0.0356
7 -6 0.0233
8 -5 0.0110
9 -4 -0.0013
10 -3 -0.0136
etc.
4) Yes, you can keep the trend model and replace the year covariates with weather factors. That's what I'd do.
Best,
Jeff
Jeffrey A. Hostetler, Ph.D.
Postdoctoral Fellow
Migratory Bird Center
Smithsonian Conservation Biology Institute
National Zoological Park
MRC 5503, Washington, DC 20013-7012
Dan Linden
As for the 2 seasons, you could expand your primary periods from 25 to 50. Then you could have year as a covariate that gets repeated for each season (e.g. 1, 1, 2, 2, 3, 3, etc.) if you thought there were yearly effects that impact both seasons equally. You then might also want season as a covariate for omega and gamma, depending on the ecology of your system. As you might imagine, things can get complicated pretty fast. But if you have good data (which it sounds like you do) and solid hypotheses, you could learn quite a bit from this work. I'd recommend getting a good book on GLMs in R before going too crazy with unmarked.
fle...@ag.arizona.edu
Dear Jeff, Dan, et al.
Thanks for the feedback gentlemen, it’s very helpful! I should have known you have to exp^ the coefficient to compute lambda given it’s a Poisson process model.
Hoping you can help with an additional question on interpretation, and clarify the data structure for the 2-season model. As suggested by Jeff, I replaced the year effect by fitting various (unstandardized) weather factors in the process model to assess hypothesized weather effects on temporal variation in population growth rate. Some model selection results with the first 2 models follow:
|
Hypothesis (predicted effect) |
nPars |
deltaAIC |
|
Prey Enhancement (ws rain t-1) |
11 |
0 |
|
Prey Enhancement (ws rain t-1) + Reproductive Stress (ws mean max. temp t-1) |
12 |
1.3 |
|
Prey Enhancement (ws rain t-1) + Winter Stress (cs mean min temp t-0.5) |
12 |
1.88 |
|
Winter Stress (cs mean min temp t-0.5) |
11 |
5.79 |
|
Reproductive Stress (ws mean max. temp t-1) |
11 |
7.01 |
|
Prey Enhancement (ws rain t-1) + Predation (mean annual rain t-2,3,4) |
12 |
17.3 |
|
Predation (mean annual rain t-2,3,4) |
11 |
22.84 |
Call:
pcountOpen(lambdaformula = ~1, gammaformula = ~Pws, omegaformula = ~1, pformula = ~Temp + I(Temp^2) + TAS + I(TAS^2) + Topo + MRain, data = cnti.umf, mixture = "P", K = 70, dynamics = "trend", se = T)
Abundance:
Estimate SE z P(>|z|)
2.61 0.131 19.9 2.69e-88
Recruitment:
Estimate SE z P(>|z|)
(Intercept) 0.12153 0.042208 2.88 0.00398
Pws -0.00125 0.000407 -3.06 0.00220
Detection:
Estimate SE z P(>|z|)
(Intercept) -2.1091 0.1668 -12.64 1.23e-36
Temp -0.1106 0.0579 -1.91 5.63e-02
I(Temp^2) -0.3185 0.0308 -10.34 4.68e-25
TAS 0.3701 0.0567 6.52 6.82e-11
I(TAS^2) -0.0660 0.0262 -2.52 1.18e-02
Toporocky slope -0.9199 0.2007 -4.58 4.57e-06
Topovalley floor 0.2688 0.1525 1.76 7.80e-02
MRain -0.0766 0.0247 -3.10 1.95e-03
AIC: 7678.935
Call:
pcountOpen(lambdaformula = ~1, gammaformula = ~Pws + Tws_max, omegaformula = ~1, pformula = ~Temp + I(Temp^2) + TAS + I(TAS^2) + Topo + MRain, data = cnti.umf, mixture = "P", K = 70, dynamics = "trend", se = T)
Abundance:
Estimate SE z P(>|z|)
2.62 0.13 20.1 3.6e-90
Recruitment:
Estimate SE z P(>|z|)
(Intercept) 1.78301 0.621129 2.87 0.004097
Pws -0.00141 0.000419 -3.36 0.000775
Tws_max -0.04563 0.016915 -2.70 0.006979
Detection:
Estimate SE z P(>|z|)
(Intercept) -2.1443 0.1701 -12.60 2.00e-36
Temp -0.1102 0.0577 -1.91 5.64e-02
I(Temp^2) -0.3144 0.0307 -10.23 1.38e-24
TAS 0.3678 0.0565 6.51 7.62e-11
I(TAS^2) -0.0657 0.0262 -2.51 1.21e-02
Toporocky slope -0.9113 0.2039 -4.47 7.82e-06
Topovalley floor 0.2748 0.1549 1.77 7.61e-02
MRain -0.0727 0.0247 -2.94 3.23e-03
AIC: 7680.236
Based on Jeff’s logic, I assume the coefficient for Pws (warm-season rain in mm 1 yr prior) in the first model means only that rainfall is declining across time, and not that each additional mm of rain has a exp^-0.00125 on the % change in abundance, as in a Poisson GLM (and the same for Tws_max in the 2nd model holding Pws constant). If that’s the case, wouldn’t the presence of trends in the explanatory variables alone contribute to changes in AIC scores of various models and bias model selection? How do I calculate the effect? Is a 100 mm increase in Pws equal to the intercept + (Coef Pws*change in Pws) rain or exp^(0.1215 + (-0.00125*100)) = 0.9965 or a 0.35% decline in population growth rate, as suggested by Jeff’s math? Biologically, increasing Pws should have a positive effect on abundance, and does based on a model using the raw counts.
I’m not sure whether I have my data structured as Dan suggested, which is important because it doubles sample sizes! In this case, I want to know how overall population growth varies across time and is affected by weather factors, and so average effects between seasons are fine. I may have formatted it correctly but am getting the following error even though all names in my code match those in the data file:
Error in `[.data.frame`(CNTI, , c("T1989_1sp", "T1989_2sp", "T1989_3sp", : undefined columns selected
To format, I placed count data for each year and season side by side (my 200 columns for 25 yrs x 8 secondary surveys in spring are next to the 25 yrs x 8 block for summer) and did the same with each ObsCov. I then entered each yearlySiteCov as 2 columns x the number of years, so there is 1 for each season in each year. So for the year covariate Dan suggested I have the column names (and values): Yr1989su (=1), Yr1990su (=2), Yr1991su (=3)…., Yr1989sp (=1), Yr1990sp (=2), Yr1991s (=3)... So 2 blocks of 25 columns side by side. I assume I then create one covariate for year as: yearlySiteCovs = list(Yr=CNTI[,c("Yr1989sp","Yr1989su","Yr1990sp","Yr1990su",…. and then fit Yr in my process model. I can provide code and data if anyone would be willing to help me.
Thanks so much folks!
Aaron
Jeff
You're welcome! I'm afraid you're misinterpreting your new results somewhat, however. The model you ran before described changes in population growth rate as a trend over time because you fit year as a numerical covariate. These new models directly model population growth rate as a function of the weather covariates. So yes, an increase of pws by 100 mm should decrease r by 0.125 or multiply lambda by exp(-0.125)=0.882. The intercept gives the estimated population growth rate (on the log scale) if pws were 0. If the effect of pws matches neither your biological intuition nor your raw data patterns, I'd suggest first checking to make sure you got your weather data into unmarked correctly. There are some potentially tricky things about the way unmarked handles data that have gotten me, at least, into trouble (e.g., R by default converts vectors into matrices in column-major order, but unmarked does it in row-major order).
Best,
Jeff
<w:LsdException Locked="false" Priority="66" SemiHidden="fa...