pcount
alejandro....@gmail.com
I am in the process of analysing some long-term monitoring data from audio-visual counts of rainforest birds. My data are spatially and temporally replicated, with "points" in "sites" distributed across a broad environmental gradient (elevation, proxied here for simplicity with mean annual temperature ("MATemp")) and repeated several times a year per site for about 10 years.
my very holey count data look like this:
count.1 count.2 count.3 count.4 count.5 count.6 count.7 count.8 count.9 count.10 count.11 count.12 count.13
KUBC3 0 2 NA NA NA NA NA NA NA NA NA NA NA
KUBC4 0 NA NA NA NA NA NA NA NA NA NA NA NA
KUBC5 1 NA NA NA NA NA NA NA NA NA NA NA NA
KUBC6 0 NA NA NA NA NA NA NA NA NA NA NA NA
TU8A2 0 NA NA NA NA NA NA NA NA NA NA NA NA
AU10A5 3 4 0 0 1 4 0 1 5 0 3 NA NA
AU10A6 2 0 0 0 0 3 4 2 2 4 3 NA NA
AU10A3 0 0 1 3 2 3 0 5 1 2 2 1 1
AU10A2 0 1 0 0 2 0 2 1 0 NA NA NA NA
siteCovs look like this
MATemp annual_pptn12
KUBC3 19.4 2120
KUBC4 19.3 2144
KUBC5 19.3 2144
KUBC6 19.3 2144
KUBC2 19.3 2153
KUBA1 20.0 1922
obsCovs look like this
wind.1 wet.1 temp_anomaly.1 month.1 start2.1 wind.2 wet.2 temp_anomaly.2 month.2 start2.2 wind.3 wet.3
KUBC3 0 1 1.6 3 8.10 NA NA NA NA NA NA NA
KUBC4 0 1 -1.3 10 7.13 1 1 -0.3 3 6.40 NA NA
KUBC5 0 1 1.2 3 7.35 NA NA NA NA NA NA NA
KUBC6 2 1 0.2 3 8.30 NA NA NA NA NA NA NA
KUBC2 0 1 1.7 3 7.25 2 1 -1.3 3 6.27 NA NA
KUBA1 2 1 1.0 10 8.25 NA NA NA NA NA NA NA
At the moment I am focussing on getting some reasonable models fitted incorporating covariates of detection and abundance, so ignoring the temporal component for now. I am able to return some good fits with a very simple quadratic term for temperature, looking like this (see above, for some reason it wants to display there...).
Based on AIC, this model performs better than one without covariates, without a quadratic terms, or with alternative error distributions...
nPars AIC delta AICwt cumltvWt
lam(I(MATemp^2))p(.)NB 4 3138.53 0.00 1.0e+00 1.00
lam(MATemp)p(.)NB 4 3163.51 24.98 3.8e-06 1.00
lam(MATemp)p(.)ZIP 4 3213.17 74.65 6.2e-17 1.00
lam(I(MATemp^2))p(.)P 3 3368.83 230.30 9.8e-51 1.00
lam(MATemp)p(.)P 3 3504.11 365.58 4.1e-80 1.00
lam(.)p(.)P 2 3633.74 495.22 2.9e-108 1.00
and has reasonable bootstrap support,
Call: parboot(object = fm6, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE 3892 -1214 1524.3 0.7624
Chisq 8212 2950 1126.6 0.0297
freemanTukey 807 -70 96.5 0.7030
But.... it returns max abundance estimates higher than 30, nearly three times the maximum recorded count (~12 individuals) (when plotted as above on the original scale). From a search of these pages and others, this could result from very low detectability estimates, and indeed, I have a lot of zeros (though this model also outperformed its ZIP equivalents) but when I try to include observation covariates to absorb some of this variation, e.g. including a covariate for start time of the survey:
> (fm8 <- pcount(~start ~I(MATemp^2), spp.umf, starts = c(1,0,0,0,0), K = 100,mixture = "NB"))
I cannot get past this error message
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
I have tried rescaling this covariate, but as it is catergorical this was possibly not even appropriate. I have five other covariates, month, wind, wet, and even temperature anomaly (the deviation of the survey from the mean temp at a site), each rescaled and raw, to no avail. there is wide variation in the number of visits to my sites, ranging from only 3 visits to 18, by when I restrict occaisions to 3 max, I get the same error. I have also tried tweaking starting values, but I am not sure I know how to choose reasonable ones in this case... Is pcount even appropriate n the case of data like these?
Any thoughts would be much appreciated!
regards,
Alex
Richard Chandler
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
Alex Anderson
Thanks so much for your quick reply. So I have tried again with the following code, adding the linear term as you suggested:
> (fm8 <- pcount(~wet ~MATemp+I(MATemp^2), spp.umf, starts = c(1,0,0,0,0,0), K = 100,mixture = "NB"))
.... and without starting values (incorrect number of starting values results in a separate specific error message)
> (fm8 <- pcount(~wet ~MATemp+I(MATemp^2), spp.umf, K = 100,mixture = "NB"))
And I still get the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
Interestingly, this same error message is also returned when I try to fit a similar model using the example scripts and data from the 2012 webinar contained in the file: moduleABUND_MHB.R:
File starts:
#mhbdata <- read.csv("http://sites.google.com/site/unmarkedinfo/home/webinars/2012-january/data/wtmatrix.csv?attredirects=0&d=1")
Later....
# Exercise
# Expand the model set to include some effects on p (duration, day, quadratics?)
fm9<- pcount(~1 ~ forest+elev+I(elev^2) + length,mhb.umf) this works...
fm10<- pcount(~day ~ forest+elev+I(elev^2) + length,mhb.umf)this doesn't...
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
cheers
A
Richard Chandler
Richard Chandler
Alex Anderson
Thanks for your reply. These models now work fine with obscovs when I run them using engine = "R", though a little slow as you expected. I am still getting very high predicted abundances... I loom forward to hearing the results of your efforts to fix the NA handling in the C++ code.
Cheers
A
Kery Marc
the NegBin often fits, but can produce unrealistically high estimates of N; see paper by Johnson et al sometime back in 2009 or so. I would clearly not use it in this case. What about the ZIP ? Does this produce reasonable estimates ?
Re. the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
Totally inexplicable to me (and to Richard Chandler as well), for about 2 months, I have had the same problem when fitting Nmix models with NAs in the data set, even with the pcount example and the mallard data set. See here:
> # Real data
> data(mallard)
> mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
+ obsCovs = mallard.obs)
> (fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
Fehler in optim(starts, nll, method = method, hessian = se, ...) :
Anfangswert in 'vmmin' ist nicht endlich
Zusätzlich: Warnmeldung:
4 sites have been discarded because of missing data.
When you fill all NAs in the covariate data, the problem goes away. Very strange.
Kind regards -- Marc
______________________________________________________________
Marc Kéry
Tel. ++41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
*** Introduction to Bayesian statistical modeling: Kéry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: Kéry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
Sent: 01 August 2013 06:39
To: unma...@googlegroups.com
Subject: [unmarked] pcount
Kery Marc
that sounds like a good explanation to me for that error.
Regards -- Marc
Marc Kéry
Tel. ++41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
*** Introduction to Bayesian statistical modeling: Kéry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: Kéry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
Sent: 02 August 2013 14:56
To: Unmarked package
Subject: Re: [unmarked] pcount
Murray Efford
Would that be Joseph, Elkin, Martin & Possingham (2009) Modeling abundance using N-mixture models: the importance of considering ecological mechanisms. Ecol. Appl. 19:631-642?
It seems to fit. I'm curious how we deal convincingly with strong model-dependence in these cases. Perhaps we can rely on the accumulated wisdom of practitioners, but that is a little hard to justify to statisticians!
Murray
Sent: Thursday, 8 August 2013 10:36 p.m.
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Kery Marc
yes, that's the paper. And several people, including myself, have fitted Nmix models with Poisson or NegBin priors for abundance and got unrealistic abundance estimates from the latter, even when a traditional GOF test (e.g., based on Chisquare) indicated the model fit. (As an aside, this is a good point to remember: whether a model fits or not does not necessarily mean anything in terms of whether it is useful.)
I find this a difficult problem: to decide which mixture to adopt for N in the model. Quite often, we find that a Poisson mixture does not fit, even when we add a couple of covariates. Since traditional wisdom says we should not base our inference on a model that does not pass some GOF test, we should therefore try some other mixture distribution, e.g., the ZIP or the NegBin which are currently implenented in unmarked. Others that have been fit in the context of a Nmix model are the Poisson log-normal or a DPP (for the latter, see the 2008 Biometrics paper by Dorazio et al.). There is clearly scope for research here.
Kind regards -- Marc
Sent: 08 August 2013 14:03
Alex Anderson
Thanks so much for sharing helpful comments. Following Richards suggestions, I've used the R engine and succeeded in having pcount fit some Nmix models with obscovs to my bird survey data from Australian montane rainforests. Today, after a package update, I am even having success with the C engine (Thanks to Richard?). As Marc and others have pointed out before, it is possible to get models with a good AIC performance and GOF that nonetheless are way above the estimates one would expect from the data. I have a set of models in which a ZIP model with no detection covariates out-performs NB, but only just.
fm13~ta~T+T^2,NB 6 2910.24 0.92 2.8e-01 0.74
fm20:~Wt+ta~T+T^2,ZIP 8 2911.83 2.50 1.3e-01 0.87
fm22:~Wt+Sn+ta~T+T^2+Pptn,ZIP 8 2911.83 2.50 1.3e-01 0.99
fm15:~Sn~T+T^2,NB 6 2920.72 11.40 1.5e-03 1.00
fm5:~1~T+T^2,NB 5 2920.86 11.54 1.4e-03 1.00
fm11~Wn~T+T^2,NB 6 2921.48 12.16 1.0e-03 1.00
fm18:~St~T+T^2,ZIP 6 2921.53 12.21 1.0e-03 1.00
fm9:~Wt~T+T^2,NB 6 2921.95 12.63 8.2e-04 1.00
fm8:~1~T+T^2+Pptn,NB 6 2948.27 38.94 1.6e-09 1.00
fm14:~ta~T+T^2,ZIP 6 2961.14 51.81 2.5e-12 1.00
fm21:~Wt+Sn+ta~T+T^2+Pptn,NB 8 2962.04 52.71 1.6e-12 1.00
fm6:~1~T+T^2,ZIP 5 2969.38 60.06 4.1e-14 1.00
fm12:~Wn~T+T^2,ZIP 6 2969.71 60.39 3.5e-14 1.00
fm19:~Wt+ta~T+T^2,NB 6 2970.11 60.79 2.9e-14 1.00
fm16:~Sn~T+T^2,ZIP 6 2970.35 61.02 2.5e-14 1.00
fm17:~St~T+T^2,NB 6 2970.35 61.02 2.5e-14 1.00
fm10:~Wt~T+T^2,ZIP 6 2971.24 61.91 1.6e-14 1.00
fm3:~1~T,NB 4 3134.72 225.40 5.1e-50 1.00
fm4:~1~T,ZIP 4 3170.93 261.61 7.0e-58 1.00
fm2:~1~T,P 3 3464.14 554.82 1.5e-121 1.00
fm1:~1~1,P 2 3591.60 682.28 3.2e-149 1.00
(site covariates: "T" = mean annual temperature, "Pptn" =mean annual precipitation,
obscovs: "ta" = temp anomaly of the survey (= survey temp minus mean annual temperature),
"Wt" = survey wetness (= mostly canopy drip in the rainforest), "Sn" = Season, "Wn" = wind,
"St" = start time)
Model fits are good for both these top models:
#fm7:~1~T+T^2,ZIP
Call: parboot(object = fm7, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
Chisq 5819 2011.9 1046.9 0.0396
freemanTukey 613 52.2 43.9 0.0891
#fm13~ta~T+T^2,NB
Call: parboot(object = fm13, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
Chisq 5529 1690.8 1268.7 0.0396
freemanTukey 603 13.6 39.7 0.3465
fm20:~Wt+ta~T+T^2,ZIP
Call: parboot(object = fm20, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
Chisq 5460 1770.4 1058.4 0.0594
freemanTukey 600 19.6 42.9 0.3168
but the estimates they give vary widely! e.g., in this case, where the observed max(count) = 12 individuals, the maximum estimate from a NegBin (top figure) is around 50 (!), while that from ZIP with obscovs is around 13.8 (middle figure), and around 20 without obscovs (bottom figure). detectability is reasonably good in this species, despite the habitat context Lewins' Honeyeater, a medium sized, vocal passerine with a loud and distinctive call)
As pointed out previously by Richard, this brings me up against a current limitation in the unmarked code for predict where the prior is a ZIP function: At the moment there are no confidence intervals calculated for a ZIP function. Confidence intervals for zero-inflated functions (Neg Bin also) (and perhaps the quicker C++ code to run them!?) would be a fantastic addition to unmarked.. In the meantime, is there any example code out there for a recommended way to achieve this outside unmarked in the mean time?
Thanks again for your helpful comments.
regards
Alex
Kery Marc
Dear Alex,
re. CIs for a ZIP fit of the Nmix, you can use a predict function in Marc Mazerolle’s AICcmodavg package (pointed out to me by Andy Royle).
Regards -- Marc
Kery Marc
Question to Alex: so you no longer got the infinite initial in vmmin error when fitting the model after you updated unmarked ? I just did and I still get it ....
Richard Chandler
Alex Anderson
Thanks for the tip for the ZIP confidence intervals I am trying that now, but seem to have a slight difference in structure of newdata than it is expecting? Everything is as for predict() in unmarked, and this command works for returning the estimates without CIs:
E.psi <- predict(fm7, type="state", newdata=newData2, appendData=TRUE)
> head(E.psi)
Predicted SE lower upper MATemp annual_pptn12
1 0.3034907 0.2007602 0.08299893 1.109733 -2.233072 -1.0716072
2 0.5255922 0.3116520 0.16440900 1.680243 -2.133072 -0.9401777
3 0.8765762 0.4634873 0.31097050 2.470928 -2.033072 -0.8087481
4 1.4078869 0.6604688 0.56137297 3.530889 -1.933072 -0.6773185
5 2.1776253 0.9024051 0.96659900 4.905914 -1.833072 -0.5458889
6 3.2436641 1.1840087 1.58609281 6.633506 -1.733072 -0.4144594
but I when I try this from package "AICcmodavg":
predictSE.zip(fm7, newdata = newData2, se.fit = TRUE, parm.type = "lambda",
type = "response", print.matrix = TRUE)
I get the error:
Error in `[.data.frame`(site.cov.data, , var.names) :
undefined columns selected
Re the infinite initial vmin problem, I am no longer getting that error, but shifting from C++ to R engine as Richard suggested, seems to have been the fix. I have updated also, but this may be unrelated? Richard?
Cheers,
Alex
Kery Marc
hi Richard,
I am using R version 2.15.2 (2012-10-26) -- "Trick or Treat" and the latest version of the packages. Perhaps I should upgrade R then ?
Regards --- Marc
Kery Marc
Dear Alex,
re. AICcmodavg: I haven’t tried out yet ...
Regards --- Marc
Alex Anderson
Until recently I have been getting good results fitting models in pcount open to bird survey data.� I then relaxed some of my data filtering to allow more variation in survey conditions, thinking that this may help to the improve fit of the detectability component of my models.� Until this point I was getting good quadratic fits to my data (see plots below), (and able to predict from ZIP fits using AICcmodavg).� Now I seem to be getting very flat or sometime even concave fits to the same data!
I have gone back to my original data now to test, and the results are similar.Adding many observation covariates (overfitting essentially) makes little difference, as does changing the family (P, NB) Have I missed something?� Attached is a zipped archive with script and files necessary to recreate the models and plot.� Any thoughts would be much appreciated.
best regards
Alex
On 9/08/13 1:34 PM, Richard Chandler wrote:
Hi Marc, what versions of unmarked and RcppArmadillo are you using?�
Regarding the delta method for the ZIP model, I'm afraid I don't have time to write the code right now, but it would be great if someone wanted to contribute code.�
Richard
On Fri, Aug 9, 2013 at 4:21 AM, Kery Marc <marc...@vogelwarte.ch> wrote:
Question to Alex: so you no longer got the infinite initial in vmmin error when fitting the model after you updated unmarked ? I just did and I still get it ....
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Kery Marc
Gesendet: Freitag, 9. August 2013 10:17
An: 'unma...@googlegroups.com'
Betreff: AW: [unmarked] pcount
�
Dear Alex,
�
re. CIs for a ZIP fit of the Nmix, you can use a predict function in Marc Mazerolle�s AICcmodavg package (pointed out to me by Andy Royle).
�
Regards� --� Marc
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Alex Anderson
Gesendet: Freitag, 9. August 2013 09:37
An: unma...@googlegroups.com
Betreff: Re: [unmarked] pcount
�
Dear all,
Thanks so much for sharing helpful comments.� Following Richards suggestions, I've used the R engine and succeeded in having pcount fit some Nmix models with obscovs to my bird survey data from Australian montane rainforests.� Today, after a package update, I am even having success with the C engine (Thanks to Richard?).� As Marc and others have pointed out before, it is possible to get models with a good AIC performance and GOF that nonetheless are way above the estimates one would expect from the data.� I have a set of models in which a ZIP model with no detection covariates out-performs NB, but only just.
���������������������������� nPars���� AIC� delta��� AICwt cumltvWt
fm7:~1~T+T^2+Pptn,ZIP������������ 6 2909.32�� 0.00� 4.5e-01���� 0.45
fm13~ta~T+T^2,NB����������������� 6 2910.24�� 0.92� 2.8e-01���� 0.74
fm20:~Wt+ta~T+T^2,ZIP������������ 8 2911.83�� 2.50� 1.3e-01���� 0.87
fm22:~Wt+Sn+ta~T+T^2+Pptn,ZIP���� 8 2911.83�� 2.50� 1.3e-01���� 0.99
fm15:~Sn~T+T^2,NB���������������� 6 2920.72� 11.40� 1.5e-03���� 1.00
fm5:~1~T+T^2,NB������������������ 5 2920.86� 11.54� 1.4e-03���� 1.00
fm11~Wn~T+T^2,NB����������������� 6 2921.48� 12.16� 1.0e-03���� 1.00
fm18:~St~T+T^2,ZIP��������������� 6 2921.53� 12.21� 1.0e-03���� 1.00
fm9:~Wt~T+T^2,NB����������������� 6 2921.95� 12.63� 8.2e-04���� 1.00
fm8:~1~T+T^2+Pptn,NB������������� 6 2948.27� 38.94� 1.6e-09���� 1.00
fm14:~ta~T+T^2,ZIP��������������� 6 2961.14� 51.81� 2.5e-12���� 1.00
fm21:~Wt+Sn+ta~T+T^2+Pptn,NB����� 8 2962.04� 52.71� 1.6e-12���� 1.00
fm6:~1~T+T^2,ZIP����������������� 5 2969.38� 60.06� 4.1e-14���� 1.00
fm12:~Wn~T+T^2,ZIP��������������� 6 2969.71� 60.39� 3.5e-14���� 1.00
fm19:~Wt+ta~T+T^2,NB������������� 6 2970.11� 60.79� 2.9e-14���� 1.00
fm16:~Sn~T+T^2,ZIP��������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm17:~St~T+T^2,NB���������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm10:~Wt~T+T^2,ZIP��������������� 6 2971.24� 61.91� 1.6e-14���� 1.00
fm3:~1~T,NB���������������������� 4 3134.72 225.40� 5.1e-50���� 1.00
fm4:~1~T,ZIP��������������������� 4 3170.93 261.61� 7.0e-58���� 1.00
fm2:~1~T,P����������������������� 3 3464.14 554.82 1.5e-121���� 1.00
fm1:~1~1,P����������������������� 2 3591.60 682.28 3.2e-149���� 1.00
(site covariates: "T" = mean annual temperature, "Pptn" =mean annual precipitation,
obscovs: "ta" = temp anomaly of the survey (= survey temp minus mean annual temperature),
"Wt" = survey wetness (= mostly canopy drip in the rainforest), "Sn" = Season, "Wn" = wind,
"St" = start time)
Model fits are good for both these top models:
#fm7:~1~T+T^2,ZIP
Call: parboot(object = fm7, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3329�������� -111.1����������� 620.4������ 0.5149
Chisq������� 5819�������� 2011.9���������� 1046.9������ 0.0396
freemanTukey� 613���������� 52.2������������ 43.9������ 0.0891
#fm13~ta~T+T^2,NB
Call: parboot(object = fm13, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3293�������� -529.9����������� 761.4������ 0.7426
Chisq������� 5529�������� 1690.8���������� 1268.7������ 0.0396
freemanTukey� 603���������� 13.6������������ 39.7������ 0.3465
fm20:~Wt+ta~T+T^2,ZIP
Call: parboot(object = fm20, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3279�������� -473.4����������� 744.1������ 0.7525
Chisq������� 5460�������� 1770.4���������� 1058.4������ 0.0594
freemanTukey� 600���������� 19.6������������ 42.9������ 0.3168
but the estimates they give vary widely! e.g., in this case, where the observed max(count)� = 12 individuals, the maximum estimate from a NegBin (top figure) is around 50 (!), while that from ZIP with obscovs is around 13.8 (middle figure), and around 20 without obscovs (bottom figure).� detectability is reasonably good in this species, despite the habitat context� Lewins' Honeyeater, a medium sized, vocal passerine with a loud and distinctive call)
As pointed out previously by Richard, this brings me up against a current limitation in the unmarked code for predict where the prior is a ZIP function:� At the moment there are no confidence intervals calculated for a ZIP function.� Confidence intervals for zero-inflated functions (Neg Bin also) (and perhaps the quicker C++ code to run them!?) would be a fantastic addition to unmarked.. In the meantime, is there any example code out there for a recommended way to achieve this outside unmarked in the mean time?
Thanks again for your helpful comments.
regards
Alex
On 8/08/13 10:40 PM, Kery Marc wrote:
Dear Murray,
yes, that's the paper. And several people, including myself, have fitted Nmix models with Poisson or NegBin priors for abundance and got unrealistic abundance estimates from the latter, even when a traditional GOF test (e.g., based on Chisquare) indicated the model fit. (As an aside, this is a good point to remember: whether a model fits or not does not necessarily mean anything in terms of whether it is useful.)
I find this a difficult problem: to decide which mixture to adopt for N in the model. Quite often, we find that a Poisson mixture does not fit, even when we add a couple of covariates. Since traditional wisdom says we should not base our inference on a model that does not pass some GOF test, we should therefore try some other mixture distribution, e.g., the ZIP or the NegBin which are currently implenented in unmarked. Others that have been fit in the context of a Nmix model are the Poisson log-normal or a DPP (for the latter, see the 2008 Biometrics paper by Dorazio et al.). There is clearly scope for research here.
Kind regards� --� Marc
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Murray Efford [murray...@otago.ac.nz]
Sent: 08 August 2013 14:03
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Hi Marc et al
Would that be Joseph, Elkin, Martin & Possingham (2009) Modeling abundance using N-mixture models: the importance of considering ecological mechanisms. Ecol. Appl. 19:631-642?
It seems to fit. I'm curious how we deal convincingly with strong model-dependence in these cases. Perhaps we can rely on the accumulated wisdom of practitioners, but that is a little hard to justify to statisticians!
Murray
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]
Sent: Thursday, 8 August 2013 10:36 p.m.
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Dear Alejandro,
the NegBin often fits, but can produce unrealistically high estimates of N; see paper by Johnson et al sometime back in 2009 or so. I would clearly not use it in this case. What about the ZIP ? Does this produce reasonable estimates ?
Re. the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
Totally inexplicable to me (and to Richard Chandler as well), for about 2 months, I have had the same problem when fitting Nmix models with NAs in the data set, even with the pcount example and the mallard data set. See here:
>�� # Real data
>����� data(mallard)
>����� mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
+����� obsCovs = mallard.obs)
>����� (fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
Fehler in optim(starts, nll, method = method, hessian = se, ...) :
� Anfangswert in 'vmmin' ist nicht endlich
Zus�tzlich: Warnmeldung:
4 sites have been discarded because of missing data.
When you fill all NAs in the covariate data, the problem goes away. Very strange.
Kind regards� --� Marc
______________________________________________________________�
Marc K�ry
�
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
�
*** Introduction to Bayesian statistical modeling: K�ry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: K�ry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of alejandro....@gmail.com [alejandro....@gmail.com]
Sent: 01 August 2013 06:39
To: unma...@googlegroups.com
Subject: [unmarked] pcount
Hi All,
I am in the process of analysing some long-term monitoring data from audio-visual counts of rainforest birds.� My data are spatially and temporally replicated, with "points" in "sites" distributed across a broad environmental gradient (elevation, proxied here for simplicity with mean annual temperature ("MATemp")) and repeated several times a year per site for about 10 years.�
my very holey count data look like this
:
�������� count.1 count.2 count.3 count.4 count.5 count.6 count.7 count.8 count.9 count.10 count.11 count.12 count.13
KUBC3��������� 0���� � 2 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC4��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC5��������� 1 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC6 ������ � 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
TU8A2��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
AU10A5�������� 3������ 4������ 0������ 0������ 1������ 4������ 0������ 1������ 5������� 0������� 3������ NA������ NA
AU10A6�������� 2������ 0������ 0������ 0������ 0������ 3������ 4������ 2������ 2������� 4������� 3������ NA������ NA
AU10A3�������� 0������ 0������ 1������ 3������ 2������ 3������ 0������ 5������ 1������� 2������� 2������� 1������� 1
AU10A2�������� 0������ 1������ 0������ 0������ 2������ 0������ 2������ 1������ 0������ NA������ NA������ NA������ NA
siteCovs look like this
�������� MATemp annual_pptn12
KUBC3����� 19.4��������� 2120
KUBC4����� 19.3��������� 2144
KUBC5����� 19.3��������� 2144
KUBC6����� 19.3��������� 2144
KUBC2����� 19.3��������� 2153
KUBA1����� 20.0��������� 1922
obsCovs look like this
����� wind.1 wet.1 temp_anomaly.1 month.1 start2.1 wind.2 wet.2 temp_anomaly.2 month.2 start2.2 wind.3 wet.3
KUBC3����� 0���� 1����������� 1.6������ 3���� 8.10���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC4����� 0���� 1���������� -1.3����� 10���� 7.13����� 1���� 1���������� -0.3������ 3���� 6.40���� NA��� NA
KUBC5����� 0���� 1����������� 1.2������ 3���� 7.35���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC6����� 2���� 1����������� 0.2������ 3���� 8.30���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC2����� 0���� 1����������� 1.7������ 3���� 7.25����� 2���� 1���������� -1.3������ 3���� 6.27���� NA��� NA
KUBA1����� 2���� 1����������� 1.0����� 10���� 8.25���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
At the moment I am focussing on getting some reasonable models fitted incorporating covariates of detection and abundance, so ignoring the temporal component for now.� I am able to return some good fits with a very simple quadratic term for temperature, looking like this (see above, for some reason it wants to display there...).
Based on AIC, this model performs better than one without covariates, without a quadratic terms, or with alternative error distributions...
���������������������� nPars���� AIC� delta��� AICwt cumltvWt
lam(I(MATemp^2))p(.)NB���� 4 3138.53�� 0.00� 1.0e+00���� 1.00
lam(MATemp)p(.)NB��������� 4 3163.51� 24.98� 3.8e-06���� 1.00
lam(MATemp)p(.)ZIP�������� 4 3213.17� 74.65� 6.2e-17���� 1.00
lam(I(MATemp^2))p(.)P����� 3 3368.83 230.30� 9.8e-51���� 1.00
lam(MATemp)p(.)P���������� 3 3504.11 365.58� 4.1e-80���� 1.00
lam(.)p(.)P��������������� 2 3633.74 495.22 2.9e-108���� 1.00
and has reasonable bootstrap support,
Call: parboot(object = fm6, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3892��������� -1214���������� 1524.3������ 0.7624
Chisq������� 8212���������� 2950���������� 1126.6������ 0.0297
freemanTukey� 807����������� -70������������ 96.5������ 0.7030
�But.... it returns max abundance estimates higher than 30, nearly three times the maximum recorded count (~12 individuals) (when plotted as above on the original scale).� From a search of these pages and others, this could result from very low detectability estimates, and indeed, I have a lot of zeros (though this model also outperformed its ZIP equivalents) but when I try to include observation covariates to absorb some of this variation, e.g. including a covariate for start time of the survey:
> (fm8 <- pcount(~start ~I(MATemp^2), spp.umf, starts = c(1,0,0,0,0), K = 100,mixture = "NB"))
I cannot get past this error message
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
I have tried rescaling this covariate, but as it is catergorical this was possibly not even appropriate.� I have five other covariates, month, wind, wet, and even temperature anomaly (the deviation of the survey from the mean temp at a site), each rescaled and raw, to no avail.� there is wide variation in the number of visits to my sites, ranging from only 3 visits to 18, by when I restrict occaisions to 3 max, I get the same error.� I have also tried tweaking starting values, but I am not sure I know how to choose reasonable ones in this case...� Is pcount even appropriate n the case of data like these?
Any thoughts would be much appreciated!
regards,
Alex
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
��
--
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.
For more options, visit https://groups.google.com/groups/opt_out.--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
Richard ChandlerUniversity of GeorgiaWarnell School of Forestry and Natural Resources
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
Richard Chandler
Hi All,
Until recently I have been getting good results fitting models in pcount open to bird survey data. I then relaxed some of my data filtering to allow more variation in survey conditions, thinking that this may help to the improve fit of the detectability component of my models. Until this point I was getting good quadratic fits to my data (see plots below), (and able to predict from ZIP fits using AICcmodavg). Now I seem to be getting very flat or sometime even concave fits to the same data!
I have gone back to my original data now to test, and the results are similar.Adding many observation covariates (overfitting essentially) makes little difference, as does changing the family (P, NB) Have I missed something? Attached is a zipped archive with script and files necessary to recreate the models and plot. Any thoughts would be much appreciated.
best regards
Alex
On 9/08/13 1:34 PM, Richard Chandler wrote:
Hi Marc, what versions of unmarked and RcppArmadillo are you using?
Regarding the delta method for the ZIP model, I'm afraid I don't have time to write the code right now, but it would be great if someone wanted to contribute code.
Richard
On Fri, Aug 9, 2013 at 4:21 AM, Kery Marc <marc...@vogelwarte.ch> wrote:
Question to Alex: so you no longer got the infinite initial in vmmin error when fitting the model after you updated unmarked ? I just did and I still get it ....
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Kery Marc
Gesendet: Freitag, 9. August 2013 10:17
An: 'unma...@googlegroups.com'
Betreff: AW: [unmarked] pcount
Dear Alex,
re. CIs for a ZIP fit of the Nmix, you can use a predict function in Marc Mazerolle’s AICcmodavg package (pointed out to me by Andy Royle).
Regards -- Marc
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Alex Anderson
Gesendet: Freitag, 9. August 2013 09:37
An: unma...@googlegroups.com
Betreff: Re: [unmarked] pcount
Dear all,
Thanks so much for sharing helpful comments. Following Richards suggestions, I've used the R engine and succeeded in having pcount fit some Nmix models with obscovs to my bird survey data from Australian montane rainforests. Today, after a package update, I am even having success with the C engine (Thanks to Richard?). As Marc and others have pointed out before, it is possible to get models with a good AIC performance and GOF that nonetheless are way above the estimates one would expect from the data. I have a set of models in which a ZIP model with no detection covariates out-performs NB, but only just.
nPars AIC delta AICwt cumltvWt
fm7:~1~T+T^2+Pptn,ZIP 6 2909.32 0.00 4.5e-01 0.45
fm13~ta~T+T^2,NB 6 2910.24 0.92 2.8e-01 0.74
fm20:~Wt+ta~T+T^2,ZIP 8 2911.83 2.50 1.3e-01 0.87
fm22:~Wt+Sn+ta~T+T^2+Pptn,ZIP 8 2911.83 2.50 1.3e-01 0.99
fm15:~Sn~T+T^2,NB 6 2920.72 11.40 1.5e-03 1.00
fm5:~1~T+T^2,NB 5 2920.86 11.54 1.4e-03 1.00
fm11~Wn~T+T^2,NB 6 2921.48 12.16 1.0e-03 1.00
fm18:~St~T+T^2,ZIP 6 2921.53 12.21 1.0e-03 1.00
fm9:~Wt~T+T^2,NB 6 2921.95 12.63 8.2e-04 1.00
fm8:~1~T+T^2+Pptn,NB 6 2948.27 38.94 1.6e-09 1.00
fm14:~ta~T+T^2,ZIP 6 2961.14 51.81 2.5e-12 1.00
fm21:~Wt+Sn+ta~T+T^2+Pptn,NB 8 2962.04 52.71 1.6e-12 1.00
fm6:~1~T+T^2,ZIP 5 2969.38 60.06 4.1e-14 1.00
fm12:~Wn~T+T^2,ZIP 6 2969.71 60.39 3.5e-14 1.00
fm19:~Wt+ta~T+T^2,NB 6 2970.11 60.79 2.9e-14 1.00
fm16:~Sn~T+T^2,ZIP 6 2970.35 61.02 2.5e-14 1.00
fm17:~St~T+T^2,NB 6 2970.35 61.02 2.5e-14 1.00
fm10:~Wt~T+T^2,ZIP 6 2971.24 61.91 1.6e-14 1.00
fm3:~1~T,NB 4 3134.72 225.40 5.1e-50 1.00
fm4:~1~T,ZIP 4 3170.93 261.61 7.0e-58 1.00
fm2:~1~T,P 3 3464.14 554.82 1.5e-121 1.00
fm1:~1~1,P 2 3591.60 682.28 3.2e-149 1.00
(site covariates: "T" = mean annual temperature, "Pptn" =mean annual precipitation,
obscovs: "ta" = temp anomaly of the survey (= survey temp minus mean annual temperature),
"Wt" = survey wetness (= mostly canopy drip in the rainforest), "Sn" = Season, "Wn" = wind,
"St" = start time)
Model fits are good for both these top models:
#fm7:~1~T+T^2,ZIP
Call: parboot(object = fm7, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE 3329 -111.1 620.4 0.5149
Chisq 5819 2011.9 1046.9 0.0396
freemanTukey 613 52.2 43.9 0.0891
#fm13~ta~T+T^2,NB
Call: parboot(object = fm13, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE 3293 -529.9 761.4 0.7426
Chisq 5529 1690.8 1268.7 0.0396
freemanTukey 603 13.6 39.7 0.3465
fm20:~Wt+ta~T+T^2,ZIP
Call: parboot(object = fm20, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE 3279 -473.4 744.1 0.7525
Chisq 5460 1770.4 1058.4 0.0594
freemanTukey 600 19.6 42.9 0.3168
but the estimates they give vary widely! e.g., in this case, where the observed max(count) = 12 individuals, the maximum estimate from a NegBin (top figure) is around 50 (!), while that from ZIP with obscovs is around 13.8 (middle figure), and around 20 without obscovs (bottom figure). detectability is reasonably good in this species, despite the habitat context Lewins' Honeyeater, a medium sized, vocal passerine with a loud and distinctive call)
As pointed out previously by Richard, this brings me up against a current limitation in the unmarked code for predict where the prior is a ZIP function: At the moment there are no confidence intervals calculated for a ZIP function. Confidence intervals for zero-inflated functions (Neg Bin also) (and perhaps the quicker C++ code to run them!?) would be a fantastic addition to unmarked.. In the meantime, is there any example code out there for a recommended way to achieve this outside unmarked in the mean time?
Thanks again for your helpful comments.
regards
Alex
On 8/08/13 10:40 PM, Kery Marc wrote:
Dear Murray,
yes, that's the paper. And several people, including myself, have fitted Nmix models with Poisson or NegBin priors for abundance and got unrealistic abundance estimates from the latter, even when a traditional GOF test (e.g., based on Chisquare) indicated the model fit. (As an aside, this is a good point to remember: whether a model fits or not does not necessarily mean anything in terms of whether it is useful.)
I find this a difficult problem: to decide which mixture to adopt for N in the model. Quite often, we find that a Poisson mixture does not fit, even when we add a couple of covariates. Since traditional wisdom says we should not base our inference on a model that does not pass some GOF test, we should therefore try some other mixture distribution, e.g., the ZIP or the NegBin which are currently implenented in unmarked. Others that have been fit in the context of a Nmix model are the Poisson log-normal or a DPP (for the latter, see the 2008 Biometrics paper by Dorazio et al.). There is clearly scope for research here.
Kind regards -- Marc
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Murray Efford [murray...@otago.ac.nz]
Sent: 08 August 2013 14:03
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Hi Marc et al
Would that be Joseph, Elkin, Martin & Possingham (2009) Modeling abundance using N-mixture models: the importance of considering ecological mechanisms. Ecol. Appl. 19:631-642?
It seems to fit. I'm curious how we deal convincingly with strong model-dependence in these cases. Perhaps we can rely on the accumulated wisdom of practitioners, but that is a little hard to justify to statisticians!
Murray
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]
Sent: Thursday, 8 August 2013 10:36 p.m.
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Dear Alejandro,
the NegBin often fits, but can produce unrealistically high estimates of N; see paper by Johnson et al sometime back in 2009 or so. I would clearly not use it in this case. What about the ZIP ? Does this produce reasonable estimates ?
Re. the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
Totally inexplicable to me (and to Richard Chandler as well), for about 2 months, I have had the same problem when fitting Nmix models with NAs in the data set, even with the pcount example and the mallard data set. See here:
> # Real data
> data(mallard)
> mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
+ obsCovs = mallard.obs)
> (fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
Fehler in optim(starts, nll, method = method, hessian = se, ...) :
Anfangswert in 'vmmin' ist nicht endlich
Zusätzlich: Warnmeldung:
4 sites have been discarded because of missing data.
When you fill all NAs in the covariate data, the problem goes away. Very strange.
Kind regards -- Marc
______________________________________________________________
Marc Kéry
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
*** Introduction to Bayesian statistical modeling: Kéry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: Kéry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of alejandro....@gmail.com [alejandro....@gmail.com]
Sent: 01 August 2013 06:39
To: unma...@googlegroups.com
Subject: [unmarked] pcount
Hi All,
I am in the process of analysing some long-term monitoring data from audio-visual counts of rainforest birds. My data are spatially and temporally replicated, with "points" in "sites" distributed across a broad environmental gradient (elevation, proxied here for simplicity with mean annual temperature ("MATemp")) and repeated several times a year per site for about 10 years.
my very holey count data look like this
:
count.1 count.2 count.3 count.4 count.5 count.6 count.7 count.8 count.9 count.10 count.11 count.12 count.13
KUBC3 0 2 NA NA NA NA NA NA NA NA NA NA NA
KUBC4 0 NA NA NA NA NA NA NA NA NA NA NA NA
KUBC5 1 NA NA NA NA NA NA NA NA NA NA NA NA
KUBC6 0 NA NA NA NA NA NA NA NA NA NA NA NA
TU8A2 0 NA NA NA NA NA NA NA NA NA NA NA NA
AU10A5 3 4 0 0 1 4 0 1 5 0 3 NA NA
AU10A6 2 0 0 0 0 3 4 2 2 4 3 NA NA
AU10A3 0 0 1 3 2 3 0 5 1 2 2 1 1
AU10A2 0 1 0 0 2 0 2 1 0 NA NA NA NA
siteCovs look like this
MATemp annual_pptn12
KUBC3 19.4 2120
KUBC4 19.3 2144
KUBC5 19.3 2144
KUBC6 19.3 2144
KUBC2 19.3 2153
KUBA1 20.0 1922
obsCovs look like this
wind.1 wet.1 temp_anomaly.1 month.1 start2.1 wind.2 wet.2 temp_anomaly.2 month.2 start2.2 wind.3 wet.3
KUBC3 0 1 1.6 3 8.10 NA NA NA NA NA NA NA
KUBC4 0 1 -1.3 10 7.13 1 1 -0.3 3 6.40 NA NA
KUBC5 0 1 1.2 3 7.35 NA NA NA NA NA NA NA
KUBC6 2 1 0.2 3 8.30 NA NA NA NA NA NA NA
KUBC2 0 1 1.7 3 7.25 2 1 -1.3 3 6.27 NA NA
KUBA1 2 1 1.0 10 8.25 NA NA NA NA NA NA NA
At the moment I am focussing on getting some reasonable models fitted incorporating covariates of detection and abundance, so ignoring the temporal component for now. I am able to return some good fits with a very simple quadratic term for temperature, looking like this (see above, for some reason it wants to display there...).
Based on AIC, this model performs better than one without covariates, without a quadratic terms, or with alternative error distributions...
nPars AIC delta AICwt cumltvWt
lam(I(MATemp^2))p(.)NB 4 3138.53 0.00 1.0e+00 1.00
lam(MATemp)p(.)NB 4 3163.51 24.98 3.8e-06 1.00
lam(MATemp)p(.)ZIP 4 3213.17 74.65 6.2e-17 1.00
lam(I(MATemp^2))p(.)P 3 3368.83 230.30 9.8e-51 1.00
lam(MATemp)p(.)P 3 3504.11 365.58 4.1e-80 1.00
lam(.)p(.)P 2 3633.74 495.22 2.9e-108 1.00
and has reasonable bootstrap support,
Call: parboot(object = fm6, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE 3892 -1214 1524.3 0.7624
Chisq 8212 2950 1126.6 0.0297
freemanTukey 807 -70 96.5 0.7030
But.... it returns max abundance estimates higher than 30, nearly three times the maximum recorded count (~12 individuals) (when plotted as above on the original scale). From a search of these pages and others, this could result from very low detectability estimates, and indeed, I have a lot of zeros (though this model also outperformed its ZIP equivalents) but when I try to include observation covariates to absorb some of this variation, e.g. including a covariate for start time of the survey:
> (fm8 <- pcount(~start ~I(MATemp^2), spp.umf, starts = c(1,0,0,0,0), K = 100,mixture = "NB"))
I cannot get past this error message
Error in optim(starts, nll, method = method, hessian = se, ...) :
initial value in 'vmmin' is not finite
I have tried rescaling this covariate, but as it is catergorical this was possibly not even appropriate. I have five other covariates, month, wind, wet, and even temperature anomaly (the deviation of the survey from the mean temp at a site), each rescaled and raw, to no avail. there is wide variation in the number of visits to my sites, ranging from only 3 visits to 18, by when I restrict occaisions to 3 max, I get the same error. I have also tried tweaking starting values, but I am not sure I know how to choose reasonable ones in this case... Is pcount even appropriate n the case of data like these?
Any thoughts would be much appreciated!
regards,
Alex
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
Alex Anderson
Thanks for your quick reply.� I understand that the fit of curves, especially to� the maximum count, will not be exact unless p=1, but as you can see from the plots in me previous posts, the curves I am getting now are much flatter than before, while detectability should be similar.� I wish I could give you more detail but all I have to go on is that I have tested both my new data and the old, and cannot seem to recreate the much tighter gaussian curves I was getting before.� I cannot rule out some error in my coding, but I am suspicious of the failure to to fit at the extremes of the distribution, where there are plenty of zeroes. Any tips for how to proceed would be much appreciated.
Cheers
A
On 9/09/13 12:08 PM, Richard Chandler wrote:
Hi Alex,
I don't have time to dig through your script right now, but you should be aware that the predicted abundance curve should not "fit" the observed data points unless p=1. However, if you think pcount has changed in some way in recent versions, I would appreciate it if you could send me some evidence of that. The function gets tested before each release and I haven't noticed any problems.
Richard
On Sun, Sep 8, 2013 at 9:09 AM, Alex Anderson <alejandro....@gmail.com> wrote:
Hi All,
Until recently I have been getting good results fitting models in pcount open to bird survey data.� I then relaxed some of my data filtering to allow more variation in survey conditions, thinking that this may help to the improve fit of the detectability component of my models.� Until this point I was getting good quadratic fits to my data (see plots below), (and able to predict from ZIP fits using AICcmodavg).� Now I seem to be getting very flat or sometime even concave fits to the same data!
I have gone back to my original data now to test, and the results are similar.Adding many observation covariates (overfitting essentially) makes little difference, as does changing the family (P, NB) Have I missed something?� Attached is a zipped archive with script and files necessary to recreate the models and plot.� Any thoughts would be much appreciated.
best regards
Alex
On 9/08/13 1:34 PM, Richard Chandler wrote:
Hi Marc, what versions of unmarked and RcppArmadillo are you using?�
Regarding the delta method for the ZIP model, I'm afraid I don't have time to write the code right now, but it would be great if someone wanted to contribute code.�
Richard
On Fri, Aug 9, 2013 at 4:21 AM, Kery Marc <marc...@vogelwarte.ch> wrote:
Question to Alex: so you no longer got the infinite initial in vmmin error when fitting the model after you updated unmarked ? I just did and I still get it ....
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Kery Marc
Gesendet: Freitag, 9. August 2013 10:17
An: 'unma...@googlegroups.com'
Betreff: AW: [unmarked] pcount
�
Dear Alex,
�
re. CIs for a ZIP fit of the Nmix, you can use a predict function in Marc Mazerolle�s AICcmodavg package (pointed out to me by Andy Royle).
�
Regards� --� Marc
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Alex Anderson
Gesendet: Freitag, 9. August 2013 09:37
An: unma...@googlegroups.com
Betreff: Re: [unmarked] pcount
�
Dear all,
Thanks so much for sharing helpful comments.� Following Richards suggestions, I've used the R engine and succeeded in having pcount fit some Nmix models with obscovs to my bird survey data from Australian montane rainforests.� Today, after a package update, I am even having success with the C engine (Thanks to Richard?).� As Marc and others have pointed out before, it is possible to get models with a good AIC performance and GOF that nonetheless are way above the estimates one would expect from the data.� I have a set of models in which a ZIP model with no detection covariates out-performs NB, but only just.
���������������������������� nPars���� AIC� delta��� AICwt cumltvWt
fm7:~1~T+T^2+Pptn,ZIP������������ 6 2909.32�� 0.00� 4.5e-01���� 0.45
fm13~ta~T+T^2,NB����������������� 6 2910.24�� 0.92� 2.8e-01���� 0.74
fm20:~Wt+ta~T+T^2,ZIP������������ 8 2911.83�� 2.50� 1.3e-01���� 0.87
fm22:~Wt+Sn+ta~T+T^2+Pptn,ZIP���� 8 2911.83�� 2.50� 1.3e-01���� 0.99
fm15:~Sn~T+T^2,NB���������������� 6 2920.72� 11.40� 1.5e-03���� 1.00
fm5:~1~T+T^2,NB������������������ 5 2920.86� 11.54� 1.4e-03���� 1.00
fm11~Wn~T+T^2,NB����������������� 6 2921.48� 12.16� 1.0e-03���� 1.00
fm18:~St~T+T^2,ZIP��������������� 6 2921.53� 12.21� 1.0e-03���� 1.00
fm9:~Wt~T+T^2,NB����������������� 6 2921.95� 12.63� 8.2e-04���� 1.00
fm8:~1~T+T^2+Pptn,NB������������� 6 2948.27� 38.94� 1.6e-09���� 1.00
fm14:~ta~T+T^2,ZIP��������������� 6 2961.14� 51.81� 2.5e-12���� 1.00
fm21:~Wt+Sn+ta~T+T^2+Pptn,NB����� 8 2962.04� 52.71� 1.6e-12���� 1.00
fm6:~1~T+T^2,ZIP����������������� 5 2969.38� 60.06� 4.1e-14���� 1.00
fm12:~Wn~T+T^2,ZIP��������������� 6 2969.71� 60.39� 3.5e-14���� 1.00
fm19:~Wt+ta~T+T^2,NB������������� 6 2970.11� 60.79� 2.9e-14���� 1.00
fm16:~Sn~T+T^2,ZIP��������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm17:~St~T+T^2,NB���������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm10:~Wt~T+T^2,ZIP��������������� 6 2971.24� 61.91� 1.6e-14���� 1.00
fm3:~1~T,NB���������������������� 4 3134.72 225.40� 5.1e-50���� 1.00
fm4:~1~T,ZIP��������������������� 4 3170.93 261.61� 7.0e-58���� 1.00
fm2:~1~T,P����������������������� 3 3464.14 554.82 1.5e-121���� 1.00
fm1:~1~1,P����������������������� 2 3591.60 682.28 3.2e-149���� 1.00
(site covariates: "T" = mean annual temperature, "Pptn" =mean annual precipitation,
obscovs: "ta" = temp anomaly of the survey (= survey temp minus mean annual temperature),
"Wt" = survey wetness (= mostly canopy drip in the rainforest), "Sn" = Season, "Wn" = wind,
"St" = start time)
Model fits are good for both these top models:
#fm7:~1~T+T^2,ZIP
Call: parboot(object = fm7, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3329�������� -111.1����������� 620.4������ 0.5149
Chisq������� 5819�������� 2011.9���������� 1046.9������ 0.0396
freemanTukey� 613���������� 52.2������������ 43.9������ 0.0891
#fm13~ta~T+T^2,NB
Call: parboot(object = fm13, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3293�������� -529.9����������� 761.4������ 0.7426
Chisq������� 5529�������� 1690.8���������� 1268.7������ 0.0396
freemanTukey� 603���������� 13.6������������ 39.7������ 0.3465
fm20:~Wt+ta~T+T^2,ZIP
Call: parboot(object = fm20, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3279�������� -473.4����������� 744.1������ 0.7525
Chisq������� 5460�������� 1770.4���������� 1058.4������ 0.0594
freemanTukey� 600���������� 19.6������������ 42.9������ 0.3168
but the estimates they give vary widely! e.g., in this case, where the observed max(count)� = 12 individuals, the maximum estimate from a NegBin (top figure) is around 50 (!), while that from ZIP with obscovs is around 13.8 (middle figure), and around 20 without obscovs (bottom figure).� detectability is reasonably good in this species, despite the habitat context� Lewins' Honeyeater, a medium sized, vocal passerine with a loud and distinctive call)
As pointed out previously by Richard, this brings me up against a current limitation in the unmarked code for predict where the prior is a ZIP function:� At the moment there are no confidence intervals calculated for a ZIP function.� Confidence intervals for zero-inflated functions (Neg Bin also) (and perhaps the quicker C++ code to run them!?) would be a fantastic addition to unmarked.. In the meantime, is there any example code out there for a recommended way to achieve this outside unmarked in the mean time?
Thanks again for your helpful comments.
regards
Alex
On 8/08/13 10:40 PM, Kery Marc wrote:
Dear Murray,
yes, that's the paper. And several people, including myself, have fitted Nmix models with Poisson or NegBin priors for abundance and got unrealistic abundance estimates from the latter, even when a traditional GOF test (e.g., based on Chisquare) indicated the model fit. (As an aside, this is a good point to remember: whether a model fits or not does not necessarily mean anything in terms of whether it is useful.)
I find this a difficult problem: to decide which mixture to adopt for N in the model. Quite often, we find that a Poisson mixture does not fit, even when we add a couple of covariates. Since traditional wisdom says we should not base our inference on a model that does not pass some GOF test, we should therefore try some other mixture distribution, e.g., the ZIP or the NegBin which are currently implenented in unmarked. Others that have been fit in the context of a Nmix model are the Poisson log-normal or a DPP (for the latter, see the 2008 Biometrics paper by Dorazio et al.). There is clearly scope for research here.
Kind regards� --� Marc
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Murray Efford [murray...@otago.ac.nz]
Sent: 08 August 2013 14:03
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Hi Marc et al
Would that be Joseph, Elkin, Martin & Possingham (2009) Modeling abundance using N-mixture models: the importance of considering ecological mechanisms. Ecol. Appl. 19:631-642?
It seems to fit. I'm curious how we deal convincingly with strong model-dependence in these cases. Perhaps we can rely on the accumulated wisdom of practitioners, but that is a little hard to justify to statisticians!
Murray
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]
Sent: Thursday, 8 August 2013 10:36 p.m.
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Dear Alejandro,
the NegBin often fits, but can produce unrealistically high estimates of N; see paper by Johnson et al sometime back in 2009 or so. I would clearly not use it in this case. What about the ZIP ? Does this produce reasonable estimates ?
Re. the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
Totally inexplicable to me (and to Richard Chandler as well), for about 2 months, I have had the same problem when fitting Nmix models with NAs in the data set, even with the pcount example and the mallard data set. See here:
>�� # Real data
>����� data(mallard)
>����� mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
+����� obsCovs = mallard.obs)
>����� (fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
Fehler in optim(starts, nll, method = method, hessian = se, ...) :
� Anfangswert in 'vmmin' ist nicht endlich
Zus�tzlich: Warnmeldung:
4 sites have been discarded because of missing data.
When you fill all NAs in the covariate data, the problem goes away. Very strange.
Kind regards� --� Marc
______________________________________________________________
�
Marc K�ry
�
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
�
*** Introduction to Bayesian statistical modeling: K�ry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: K�ry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of alejandro....@gmail.com [alejandro....@gmail.com]
Sent: 01 August 2013 06:39
To: unma...@googlegroups.com
Subject: [unmarked] pcount
Hi All,
I am in the process of analysing some long-term monitoring data from audio-visual counts of rainforest birds.� My data are spatially and temporally replicated, with "points" in "sites" distributed across a broad environmental gradient (elevation, proxied here for simplicity with mean annual temperature ("MATemp")) and repeated several times a year per site for about 10 years.�
my very holey count data look like this
:
�������� count.1 count.2 count.3 count.4 count.5 count.6 count.7 count.8 count.9 count.10 count.11 count.12 count.13
KUBC3��������� 0���� � 2 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC4��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC5��������� 1 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC6 ������ � 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
TU8A2��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
AU10A5�������� 3������ 4������ 0������ 0������ 1������ 4������ 0������ 1������ 5������� 0������� 3������ NA������ NA
AU10A6�������� 2������ 0������ 0������ 0������ 0������ 3������ 4������ 2������ 2������� 4������� 3������ NA������ NA
AU10A3�������� 0������ 0������ 1������ 3������ 2������ 3������ 0������ 5������ 1������� 2������� 2������� 1������� 1
AU10A2�������� 0������ 1������ 0������ 0������ 2������ 0������ 2������ 1������ 0������ NA������ NA������ NA������ NA
siteCovs look like this
�������� MATemp annual_pptn12
KUBC3����� 19.4��������� 2120
KUBC4����� 19.3��������� 2144
KUBC5����� 19.3��������� 2144
KUBC6����� 19.3��������� 2144
KUBC2����� 19.3��������� 2153
KUBA1����� 20.0��������� 1922
obsCovs look like this
����� wind.1 wet.1 temp_anomaly.1 month.1 start2.1 wind.2 wet.2 temp_anomaly.2 month.2 start2.2 wind.3 wet.3
KUBC3����� 0���� 1����������� 1.6������ 3���� 8.10���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC4����� 0���� 1���������� -1.3����� 10���� 7.13����� 1���� 1���������� -0.3������ 3���� 6.40���� NA��� NA
KUBC5����� 0���� 1����������� 1.2������ 3���� 7.35���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC6����� 2���� 1����������� 0.2������ 3���� 8.30���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC2����� 0���� 1����������� 1.7������ 3���� 7.25����� 2���� 1���������� -1.3������ 3���� 6.27���� NA��� NA
KUBA1����� 2���� 1����������� 1.0����� 10���� 8.25���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
At the moment I am focussing on getting some reasonable models fitted incorporating covariates of detection and abundance, so ignoring the temporal component for now.� I am able to return some good fits with a very simple quadratic term for temperature, looking like this (see above, for some reason it wants to display there...).
Based on AIC, this model performs better than one without covariates, without a quadratic terms, or with alternative error distributions...
���������������������� nPars���� AIC� delta��� AICwt cumltvWt
lam(I(MATemp^2))p(.)NB���� 4 3138.53�� 0.00� 1.0e+00���� 1.00
lam(MATemp)p(.)NB��������� 4 3163.51� 24.98� 3.8e-06���� 1.00
lam(MATemp)p(.)ZIP�������� 4 3213.17� 74.65� 6.2e-17���� 1.00
lam(I(MATemp^2))p(.)P����� 3 3368.83 230.30� 9.8e-51���� 1.00
lam(MATemp)p(.)P���������� 3 3504.11 365.58� 4.1e-80���� 1.00
lam(.)p(.)P��������������� 2 3633.74 495.22 2.9e-108���� 1.00
and has reasonable bootstrap support,
Call: parboot(object = fm6, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3892��������� -1214���������� 1524.3������ 0.7624
Chisq������� 8212���������� 2950���������� 1126.6������ 0.0297
freemanTukey� 807����������� -70������������ 96.5������ 0.7030
�But.... it returns max abundance estimates higher than 30, nearly three times the maximum recorded count (~12 individuals) (when plotted as above on the original scale).� From a search of these pages and others, this could result from very low detectability estimates, and indeed, I have a lot of zeros (though this model also outperformed its ZIP equivalents) but when I try to include observation covariates to absorb some of this variation, e.g. including a covariate for start time of the survey:
> (fm8 <- pcount(~start ~I(MATemp^2), spp.umf, starts = c(1,0,0,0,0), K = 100,mixture = "NB"))
I cannot get past this error message
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
I have tried rescaling this covariate, but as it is catergorical this was possibly not even appropriate.� I have five other covariates, month, wind, wet, and even temperature anomaly (the deviation of the survey from the mean temp at a site), each rescaled and raw, to no avail.� there is wide variation in the number of visits to my sites, ranging from only 3 visits to 18, by when I restrict occaisions to 3 max, I get the same error.� I have also tried tweaking starting values, but I am not sure I know how to choose reasonable ones in this case...� Is pcount even appropriate n the case of data like these?
Any thoughts would be much appreciated!
regards,
Alex
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
��
--
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.
For more options, visit https://groups.google.com/groups/opt_out.--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
Richard ChandlerUniversity of GeorgiaWarnell School of Forestry and Natural Resources
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
Alex Anderson
Thanks for your reply Richard.� The problem I was having with strangely flat model fits was indeed due to problem with my data.� After much struggle with my data preparation code, I managed to isolate a couple of pretty simple errors, but ones that as a novice user of unmarked I found difficult to trace initially.� One of these is unique to my data, and not worth posting, but I thought to post the solution to the other in case it saves someone else the same hassle?
Basically, I had not sorted my data by site and date before creating my occasion vector, thinking that the loop below would do the trick.
for(point in unique(x$point_time_step2)){
� x$occasion[which(x$point_time_step2 == point)] = 1:length(x$occasion[which(x$point_time_step2 == point)])
}
But this created a disordered vector, and the resulting obs.covs matrix was faulty.� A simple error, but hard to find unless you examine the obs.covs matrix carefully.
I doubt its worth providing a reproducible example for such a simple mistake, but below are some data and code to explain in more detail what I did wrong:
###############################################################
My data appear something like this, with " point_time_step2" being� a vector of survey sites,� "accuMATemp". Being an example of a site covariate, "date" and "temp anomaly" being example observation covariates, and "LEWHE" an example of species count data.
> head(x,10)
�� point_time_step2 accuMATemp������ date temp_anomaly LEWHE
1��������� BK16A1_2�� 14.59313 2010-02-24���� 2.406875���� 0
2��������� BK16A1_2�� 14.59313 2009-04-19��� -2.093125���� 0
3��������� BK16A1_2�� 14.59313 2009-04-18��� -0.593125���� 0
4���������� BKAA1_1�� 14.30339 2003-10-15���� 3.196611���� 0
5���������� BKAA1_1�� 14.30339 2000-01-11���� 1.696611���� 0
6��������� BK16A1_2�� 14.59313 2009-04-20��� -2.593125���� 0
7��������� BK16A1_2�� 14.59313 2009-02-10���� 0.000000���� 0
8���������� BKAA1_1�� 14.30339 2003-10-15���� 1.696611���� 0
9��������� BK16A2_2�� 14.59313 2010-02-25���� 0.000000���� 0
10�������� BK16A2_2�� 14.59313 2010-02-25���� 4.406875���� 0
During data preparation I used the following for() loop to create an "occasion" vector to identify visits to sites (required later to create an obs covs object:
obs.x = x[,c("point_time_step2","occasion","temp_anomaly","date")]
#reshape the obs covs data into wide format with columns for visits
obs.y = reshape(obs.x, timevar = "occasion",
��������������� idvar = "point_time_step2", ids = 1:NROW(data),
��������������� times = seq_along(varying[[1]]),
��������������� drop = NULL, direction = "wide", new.row.names = NULL)
But unless I did the following:
> x = x[order(x$point_ID_2, x$date),]
��� point_time_step2 accuMATemp������ date temp_anomaly LEWHE
183�������� AU10A1_1��� 17.8926 1998-02-13��� 4.6074028���� 1
177�������� AU10A1_1��� 17.8926 1998-10-27��� 3.1074028���� 1
180�������� AU10A1_1��� 17.8926 2004-02-01��� 4.1074028���� 0
182�������� AU10A1_1��� 17.8926 2004-05-14�� -3.8925972���� 0
179�������� AU10A1_1��� 17.8926 2005-02-14��� 5.1074028���� 1
175�������� AU10A1_1��� 17.8926 2007-03-27��� 1.1074028���� 1
173�������� AU10A1_1��� 17.8926 2008-04-04�� -2.3925972���� 1
174�������� AU10A1_2��� 17.8926 2008-10-19�� -0.3925972���� 0
181�������� AU10A1_2��� 17.8926 2009-06-20�� -5.8925972���� 2
184�������� AU10A1_2��� 17.8926 2010-10-15�� -1.8925972���� 0
The occasion vector was disordered, and I think this was at least part of my problem.
The resulting fits now make much more sense!
The other major issue I had was that I created the occasion vector using only location and date, but was trying to split on time-step also in the model fitting, creating further havoc. There you go...
Cheers
A
On 9/09/13 1:08 PM, Richard Chandler wrote:
Hi Alex,
I don't have time to dig through your script right now, but you should be aware that the predicted abundance curve should not "fit" the observed data points unless p=1. However, if you think pcount has changed in some way in recent versions, I would appreciate it if you could send me some evidence of that. The function gets tested before each release and I haven't noticed any problems.
Richard
On Sun, Sep 8, 2013 at 9:09 AM, Alex Anderson <alejandro....@gmail.com> wrote:
Hi All,
Until recently I have been getting good results fitting models in pcount open to bird survey data.� I then relaxed some of my data filtering to allow more variation in survey conditions, thinking that this may help to the improve fit of the detectability component of my models.� Until this point I was getting good quadratic fits to my data (see plots below), (and able to predict from ZIP fits using AICcmodavg).� Now I seem to be getting very flat or sometime even concave fits to the same data!
I have gone back to my original data now to test, and the results are similar.Adding many observation covariates (overfitting essentially) makes little difference, as does changing the family (P, NB) Have I missed something?� Attached is a zipped archive with script and files necessary to recreate the models and plot.� Any thoughts would be much appreciated.
best regards
Alex
On 9/08/13 1:34 PM, Richard Chandler wrote:
Hi Marc, what versions of unmarked and RcppArmadillo are you using?�
Regarding the delta method for the ZIP model, I'm afraid I don't have time to write the code right now, but it would be great if someone wanted to contribute code.�
Richard
On Fri, Aug 9, 2013 at 4:21 AM, Kery Marc <marc...@vogelwarte.ch> wrote:
Question to Alex: so you no longer got the infinite initial in vmmin error when fitting the model after you updated unmarked ? I just did and I still get it ....
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Kery Marc
Gesendet: Freitag, 9. August 2013 10:17
An: 'unma...@googlegroups.com'
Betreff: AW: [unmarked] pcount
�
Dear Alex,
�
re. CIs for a ZIP fit of the Nmix, you can use a predict function in Marc Mazerolle�s AICcmodavg package (pointed out to me by Andy Royle).
�
Regards� --� Marc
�
�
�
Von: unma...@googlegroups.com [mailto:unma...@googlegroups.com] Im Auftrag von Alex Anderson
Gesendet: Freitag, 9. August 2013 09:37
An: unma...@googlegroups.com
Betreff: Re: [unmarked] pcount
�
Dear all,
Thanks so much for sharing helpful comments.� Following Richards suggestions, I've used the R engine and succeeded in having pcount fit some Nmix models with obscovs to my bird survey data from Australian montane rainforests.� Today, after a package update, I am even having success with the C engine (Thanks to Richard?).� As Marc and others have pointed out before, it is possible to get models with a good AIC performance and GOF that nonetheless are way above the estimates one would expect from the data.� I have a set of models in which a ZIP model with no detection covariates out-performs NB, but only just.
���������������������������� nPars���� AIC� delta��� AICwt cumltvWt
fm7:~1~T+T^2+Pptn,ZIP������������ 6 2909.32�� 0.00� 4.5e-01���� 0.45
fm13~ta~T+T^2,NB����������������� 6 2910.24�� 0.92� 2.8e-01���� 0.74
fm20:~Wt+ta~T+T^2,ZIP������������ 8 2911.83�� 2.50� 1.3e-01���� 0.87
fm22:~Wt+Sn+ta~T+T^2+Pptn,ZIP���� 8 2911.83�� 2.50� 1.3e-01���� 0.99
fm15:~Sn~T+T^2,NB���������������� 6 2920.72� 11.40� 1.5e-03���� 1.00
fm5:~1~T+T^2,NB������������������ 5 2920.86� 11.54� 1.4e-03���� 1.00
fm11~Wn~T+T^2,NB����������������� 6 2921.48� 12.16� 1.0e-03���� 1.00
fm18:~St~T+T^2,ZIP��������������� 6 2921.53� 12.21� 1.0e-03���� 1.00
fm9:~Wt~T+T^2,NB����������������� 6 2921.95� 12.63� 8.2e-04���� 1.00
fm8:~1~T+T^2+Pptn,NB������������� 6 2948.27� 38.94� 1.6e-09���� 1.00
fm14:~ta~T+T^2,ZIP��������������� 6 2961.14� 51.81� 2.5e-12���� 1.00
fm21:~Wt+Sn+ta~T+T^2+Pptn,NB����� 8 2962.04� 52.71� 1.6e-12���� 1.00
fm6:~1~T+T^2,ZIP����������������� 5 2969.38� 60.06� 4.1e-14���� 1.00
fm12:~Wn~T+T^2,ZIP��������������� 6 2969.71� 60.39� 3.5e-14���� 1.00
fm19:~Wt+ta~T+T^2,NB������������� 6 2970.11� 60.79� 2.9e-14���� 1.00
fm16:~Sn~T+T^2,ZIP��������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm17:~St~T+T^2,NB���������������� 6 2970.35� 61.02� 2.5e-14���� 1.00
fm10:~Wt~T+T^2,ZIP��������������� 6 2971.24� 61.91� 1.6e-14���� 1.00
fm3:~1~T,NB���������������������� 4 3134.72 225.40� 5.1e-50���� 1.00
fm4:~1~T,ZIP��������������������� 4 3170.93 261.61� 7.0e-58���� 1.00
fm2:~1~T,P����������������������� 3 3464.14 554.82 1.5e-121���� 1.00
fm1:~1~1,P����������������������� 2 3591.60 682.28 3.2e-149���� 1.00
(site covariates: "T" = mean annual temperature, "Pptn" =mean annual precipitation,
obscovs: "ta" = temp anomaly of the survey (= survey temp minus mean annual temperature),
"Wt" = survey wetness (= mostly canopy drip in the rainforest), "Sn" = Season, "Wn" = wind,
"St" = start time)
Model fits are good for both these top models:
#fm7:~1~T+T^2,ZIP
Call: parboot(object = fm7, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3329�������� -111.1����������� 620.4������ 0.5149
Chisq������� 5819�������� 2011.9���������� 1046.9������ 0.0396
freemanTukey� 613���������� 52.2������������ 43.9������ 0.0891
#fm13~ta~T+T^2,NB
Call: parboot(object = fm13, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3293�������� -529.9����������� 761.4������ 0.7426
Chisq������� 5529�������� 1690.8���������� 1268.7������ 0.0396
freemanTukey� 603���������� 13.6������������ 39.7������ 0.3465
fm20:~Wt+ta~T+T^2,ZIP
Call: parboot(object = fm20, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3279�������� -473.4����������� 744.1������ 0.7525
Chisq������� 5460�������� 1770.4���������� 1058.4������ 0.0594
freemanTukey� 600���������� 19.6������������ 42.9������ 0.3168
but the estimates they give vary widely! e.g., in this case, where the observed max(count)� = 12 individuals, the maximum estimate from a NegBin (top figure) is around 50 (!), while that from ZIP with obscovs is around 13.8 (middle figure), and around 20 without obscovs (bottom figure).� detectability is reasonably good in this species, despite the habitat context� Lewins' Honeyeater, a medium sized, vocal passerine with a loud and distinctive call)
As pointed out previously by Richard, this brings me up against a current limitation in the unmarked code for predict where the prior is a ZIP function:� At the moment there are no confidence intervals calculated for a ZIP function.� Confidence intervals for zero-inflated functions (Neg Bin also) (and perhaps the quicker C++ code to run them!?) would be a fantastic addition to unmarked.. In the meantime, is there any example code out there for a recommended way to achieve this outside unmarked in the mean time?
Thanks again for your helpful comments.
regards
Alex
On 8/08/13 10:40 PM, Kery Marc wrote:
Dear Murray,
yes, that's the paper. And several people, including myself, have fitted Nmix models with Poisson or NegBin priors for abundance and got unrealistic abundance estimates from the latter, even when a traditional GOF test (e.g., based on Chisquare) indicated the model fit. (As an aside, this is a good point to remember: whether a model fits or not does not necessarily mean anything in terms of whether it is useful.)
I find this a difficult problem: to decide which mixture to adopt for N in the model. Quite often, we find that a Poisson mixture does not fit, even when we add a couple of covariates. Since traditional wisdom says we should not base our inference on a model that does not pass some GOF test, we should therefore try some other mixture distribution, e.g., the ZIP or the NegBin which are currently implenented in unmarked. Others that have been fit in the context of a Nmix model are the Poisson log-normal or a DPP (for the latter, see the 2008 Biometrics paper by Dorazio et al.). There is clearly scope for research here.
Kind regards� --� Marc
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Murray Efford [murray...@otago.ac.nz]
Sent: 08 August 2013 14:03
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Hi Marc et al
Would that be Joseph, Elkin, Martin & Possingham (2009) Modeling abundance using N-mixture models: the importance of considering ecological mechanisms. Ecol. Appl. 19:631-642?
It seems to fit. I'm curious how we deal convincingly with strong model-dependence in these cases. Perhaps we can rely on the accumulated wisdom of practitioners, but that is a little hard to justify to statisticians!
Murray
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]
Sent: Thursday, 8 August 2013 10:36 p.m.
To: unma...@googlegroups.com
Subject: RE: [unmarked] pcount
Dear Alejandro,
the NegBin often fits, but can produce unrealistically high estimates of N; see paper by Johnson et al sometime back in 2009 or so. I would clearly not use it in this case. What about the ZIP ? Does this produce reasonable estimates ?
Re. the error message:
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
Totally inexplicable to me (and to Richard Chandler as well), for about 2 months, I have had the same problem when fitting Nmix models with NAs in the data set, even with the pcount example and the mallard data set. See here:
>�� # Real data
>����� data(mallard)
>����� mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
+����� obsCovs = mallard.obs)
>����� (fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
Fehler in optim(starts, nll, method = method, hessian = se, ...) :
� Anfangswert in 'vmmin' ist nicht endlich
Zus�tzlich: Warnmeldung:
4 sites have been discarded because of missing data.
When you fill all NAs in the covariate data, the problem goes away. Very strange.
Kind regards� --� Marc
______________________________________________________________
�
Marc K�ry
�
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________
�
*** Introduction to Bayesian statistical modeling: K�ry (2010), Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook
*** Book on Bayesian statistical modeling: K�ry & Schaub (2012), Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
*** Upcoming workshops: http://www.phidot.org/forum/viewforum.php?f=8
From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of alejandro....@gmail.com [alejandro....@gmail.com]
Sent: 01 August 2013 06:39
To: unma...@googlegroups.com
Subject: [unmarked] pcount
Hi All,
I am in the process of analysing some long-term monitoring data from audio-visual counts of rainforest birds.� My data are spatially and temporally replicated, with "points" in "sites" distributed across a broad environmental gradient (elevation, proxied here for simplicity with mean annual temperature ("MATemp")) and repeated several times a year per site for about 10 years.�
my very holey count data look like this
:
�������� count.1 count.2 count.3 count.4 count.5 count.6 count.7 count.8 count.9 count.10 count.11 count.12 count.13
KUBC3��������� 0���� � 2 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC4��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC5��������� 1 ���� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
KUBC6 ������ � 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
TU8A2��������� 0����� NA����� NA����� NA����� NA����� NA����� NA����� NA����� NA������ NA������ NA������ NA������ NA
AU10A5�������� 3������ 4������ 0������ 0������ 1������ 4������ 0������ 1������ 5������� 0������� 3������ NA������ NA
AU10A6�������� 2������ 0������ 0������ 0������ 0������ 3������ 4������ 2������ 2������� 4������� 3������ NA������ NA
AU10A3�������� 0������ 0������ 1������ 3������ 2������ 3������ 0������ 5������ 1������� 2������� 2������� 1������� 1
AU10A2�������� 0������ 1������ 0������ 0������ 2������ 0������ 2������ 1������ 0������ NA������ NA������ NA������ NA
siteCovs look like this
�������� MATemp annual_pptn12
KUBC3����� 19.4��������� 2120
KUBC4����� 19.3��������� 2144
KUBC5����� 19.3��������� 2144
KUBC6����� 19.3��������� 2144
KUBC2����� 19.3��������� 2153
KUBA1����� 20.0��������� 1922
obsCovs look like this
����� wind.1 wet.1 temp_anomaly.1 month.1 start2.1 wind.2 wet.2 temp_anomaly.2 month.2 start2.2 wind.3 wet.3
KUBC3����� 0���� 1����������� 1.6������ 3���� 8.10���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC4����� 0���� 1���������� -1.3����� 10���� 7.13����� 1���� 1���������� -0.3������ 3���� 6.40���� NA��� NA
KUBC5����� 0���� 1����������� 1.2������ 3���� 7.35���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC6����� 2���� 1����������� 0.2������ 3���� 8.30���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
KUBC2����� 0���� 1����������� 1.7������ 3���� 7.25����� 2���� 1���������� -1.3������ 3���� 6.27���� NA��� NA
KUBA1����� 2���� 1����������� 1.0����� 10���� 8.25���� NA��� NA������������ NA����� NA������ NA���� NA��� NA
At the moment I am focussing on getting some reasonable models fitted incorporating covariates of detection and abundance, so ignoring the temporal component for now.� I am able to return some good fits with a very simple quadratic term for temperature, looking like this (see above, for some reason it wants to display there...).
Based on AIC, this model performs better than one without covariates, without a quadratic terms, or with alternative error distributions...
���������������������� nPars���� AIC� delta��� AICwt cumltvWt
lam(I(MATemp^2))p(.)NB���� 4 3138.53�� 0.00� 1.0e+00���� 1.00
lam(MATemp)p(.)NB��������� 4 3163.51� 24.98� 3.8e-06���� 1.00
lam(MATemp)p(.)ZIP�������� 4 3213.17� 74.65� 6.2e-17���� 1.00
lam(I(MATemp^2))p(.)P����� 3 3368.83 230.30� 9.8e-51���� 1.00
lam(MATemp)p(.)P���������� 3 3504.11 365.58� 4.1e-80���� 1.00
lam(.)p(.)P��������������� 2 3633.74 495.22 2.9e-108���� 1.00
and has reasonable bootstrap support,
Call: parboot(object = fm6, statistic = fitstats, nsim = 100, report = 1)
Parametric Bootstrap Statistics:
�������������� t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
SSE��������� 3892��������� -1214���������� 1524.3������ 0.7624
Chisq������� 8212���������� 2950���������� 1126.6������ 0.0297
freemanTukey� 807����������� -70������������ 96.5������ 0.7030
�But.... it returns max abundance estimates higher than 30, nearly three times the maximum recorded count (~12 individuals) (when plotted as above on the original scale).� From a search of these pages and others, this could result from very low detectability estimates, and indeed, I have a lot of zeros (though this model also outperformed its ZIP equivalents) but when I try to include observation covariates to absorb some of this variation, e.g. including a covariate for start time of the survey:
> (fm8 <- pcount(~start ~I(MATemp^2), spp.umf, starts = c(1,0,0,0,0), K = 100,mixture = "NB"))
I cannot get past this error message
Error in optim(starts, nll, method = method, hessian = se, ...) :
� initial value in 'vmmin' is not finite
I have tried rescaling this covariate, but as it is catergorical this was possibly not even appropriate.� I have five other covariates, month, wind, wet, and even temperature anomaly (the deviation of the survey from the mean temp at a site), each rescaled and raw, to no avail.� there is wide variation in the number of visits to my sites, ranging from only 3 visits to 18, by when I restrict occaisions to 3 max, I get the same error.� I have also tried tweaking starting values, but I am not sure I know how to choose reasonable ones in this case...� Is pcount even appropriate n the case of data like these?
Any thoughts would be much appreciated!
regards,
Alex
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
��
--
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.
For more options, visit https://groups.google.com/groups/opt_out.--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
--
Richard ChandlerUniversity of GeorgiaWarnell School of Forestry and Natural Resources
--
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.
For more options, visit https://groups.google.com/groups/opt_out.
�
�
--
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.
For more options, visit https://groups.google.com/groups/opt_out.