colext: differences in unmarked vs. PRESENCE estimates

[Gavin M. Jones]

Hi unmarked users,

 

I'm a recent convert to unmarked from PRESENCE, and in the changeover I wanted to "check" my work to make sure I could actually do what I was meaning to do in unmarked. I originally plugged in my entire dataset and ran two simple models:

 

m1 <- colext(~1, ~1, ~1, ~1, data=umf)

m2 <- colext(~1, ~year-1, ~year-1, ~year-1, data=umf)

 

The comparison of unmarked & PRESENCE outputs is in the tables below. Even with a quick glance, there are some noted differences. AIC values are extraordinarily different (1528.036 vs. 3221.0195 in the first example), and parameter estimates are similar, but not identical.

 

The first model (m1 <- colext(~1, ~1, ~1, ~1, data=umf).

unmarked > M <- 40 > J <- 10 > T <- 20 > y <- as.matrix(mydata[,1:200]) > year <- data.frame(matrix(rep(1993:2012, each=M), M, T)) > year <- data.frame(lapply(year, as.factor)) > umf <- unmarkedMultFrame(y = y, +                          yearlySiteCovs = list(year = year), +                          numPrimary=20) > m1 <- colext(~1, ~1, ~1, ~1, data=umf) > m1   Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,     pformula = ~1, data = umf)   Initial:  Estimate  SE      z P(>|z|)      12.3 195 0.0627    0.95   Colonization:  Estimate    SE     z  P(>|z|)     -1.66 0.246 -6.72 1.77e-11   Extinction:  Estimate    SE     z  P(>|z|)     -2.35 0.189 -12.5 1.01e-35   Detection:  Estimate     SE    z  P(>|z|)      1.21 0.0792 15.3 7.81e-53   AIC: 1528.036

PRESENCE   Model(1):psi,gamma(.),eps(.),p(.)   Number of parameters           = 4 Number of function calls           = 284 -2log(likelihood)              = 3213.0195 AIC                            = 3221.0195 LikeNRSig=6 eps=0.01 ETA=1e-013   Untransformed Estimates of coefficients for covariates (Beta's) =============================================                                                       estimate    std.error A1   psi1                             :  25.064276  43543.2085       B1   gam1                             :  -0.948506    0.177658       C1   eps1                             :  -2.077670    0.134999       D1   P[1-1]                           :   1.289211    0.049367

 

 

*So, they are similar, but why the small differences in parameter estimates? Is this a user error, or am I missing something about the differences between how these two programs estimate parameters?

 

*Below, I allow ext, col, and p to vary among years. This is where the rather large differences in results are seen.

 

The second model, m2 <- colext(~1, ~year-1, ~year-1, ~year-1, data=umf)

unmarked m2 <- colext(~1, ~year-1, ~year-1, ~year-1, data=umf) > m2   Call: colext(psiformula = ~1, gammaformula = ~year - 1, epsilonformula = ~year -     1, pformula = ~year - 1, data = umf)                 Initial:  Estimate  SE   z  P(>|z|)      2.96 0.8 3.7 0.000215   Colonization:          Estimate      SE       z P(>|z|) year1993 -10.6141 151.034 -0.0703  0.9440 year1994  -5.5133     NaN     NaN     NaN year1995  -8.1765     NaN     NaN     NaN year1996  -0.0888   1.534 -0.0579  0.9538 year1997  -7.1978  46.655 -0.1543  0.8774 year1998  -7.3970  44.654 -0.1656  0.8684 year1999  -7.2636  33.214 -0.2187  0.8269 year2000  -0.8831   0.854 -1.0342  0.3010 year2001  -1.3792   1.082 -1.2753  0.2022 year2002  -0.7610   0.771 -0.9865  0.3239 year2003  -1.1307   0.924 -1.2235  0.2211 year2004  -1.5019   1.137 -1.3206  0.1867 year2005  -2.0309   1.619 -1.2541  0.2098 year2006  -1.1080   0.813 -1.3634  0.1727 year2007  -1.5113   1.082 -1.3972  0.1624 year2008  -2.9879   2.694 -1.1093  0.2673 year2009  -1.8110   0.778 -2.3270  0.0200 year2010  -9.5782  37.022 -0.2587  0.7959 year2011  -2.4318   1.184 -2.0546  0.0399   Extinction:          Estimate     SE      z P(>|z|) year1993    -9.09 35.570 -0.256 0.79820 year1994    -3.94  2.655 -1.482 0.13830 year1995    -9.85 53.596 -0.184 0.85414 year1996   -12.52 92.651 -0.135 0.89247 year1997    -3.12  1.249 -2.494 0.01262 year1998    -1.62  0.635 -2.556 0.01059 year1999    -1.40  0.648 -2.163 0.03052 year2000    -2.54  1.424 -1.786 0.07413 year2001    -1.96  0.707 -2.766 0.00567 year2002   -10.89 83.945 -0.130 0.89676 year2003    -1.01  0.484 -2.082 0.03735 year2004    -1.90  0.834 -2.277 0.02279 year2005    -2.82  2.503 -1.128 0.25938 year2006    -2.00  0.985 -2.026 0.04273 year2007    -2.27  0.926 -2.456 0.01406 year2008    -1.87  0.634 -2.944 0.00324 year2009   -10.24 54.781 -0.187 0.85178 year2010    -9.22 31.458 -0.293 0.76943 year2011    -1.30  0.556 -2.345 0.01902   Detection:          Estimate    SE     z  P(>|z|) year1993    0.979 0.442 2.214 2.69e-02 year1994    1.997 0.489 4.087 4.38e-05 year1995    1.180 0.552 2.138 3.25e-02 year1996    1.325 0.431 3.071 2.13e-03 year1997    2.420 0.394 6.136 8.46e-10 year1998    1.098 0.322 3.410 6.51e-04 year1999    1.002 0.320 3.129 1.75e-03 year2000    1.963 0.500 3.925 8.67e-05 year2001    1.444 0.362 3.993 6.52e-05 year2002    1.786 0.399 4.471 7.80e-06 year2003    1.148 0.296 3.871 1.08e-04 year2004    1.351 0.284 4.764 1.90e-06 year2005    0.962 0.380 2.534 1.13e-02 year2006    0.165 0.292 0.563 5.73e-01 year2007    1.447 0.314 4.603 4.16e-06 year2008    0.788 0.348 2.267 2.34e-02 year2009    0.695 0.285 2.436 1.48e-02 year2010    0.791 0.292 2.707 6.79e-03 year2011    0.806 0.382 2.110 3.49e-02 year2012    1.331 0.349 3.808 1.40e-04   AIC: 1558.539 Warning message: In sqrt(diag(vcov(obj))

PRESENCE Model(1):psi,gamma(year),eps(year),p(year)   Number of parameters           = 59 Number of function calls           = 5806   **** Numerical convergence may not have been reached.      Parameter estimates converged to approximately      4.05 significant digits.   -2log(likelihood)              = 3090.3353 AIC                            = 3208.3353 LikeNRSig=6 eps=0.01 ETA=1e-013   Untransformed Estimates of coefficients for covariates (Beta's)                                             estimate    std.error A1   psi1                             :  24.511994  33473.363879       B1   gam1                             :   0.198692  23512.906201 B2   gam2                             : -67.745272   10.000000 B3   gam3                             : -68.507467   10.000000 B4   gam4                             :  -0.892897    1.315030 B5   gam5                             :   0.238720    1.315971 B6   gam6                             :   0.224209    1.151735 B7   gam7                             :  -0.245802    1.028310 B8   gam8                             :  -1.501903    1.146757 B9   gam9                             :  -1.131734    0.870562 B10  gam10                            :  -0.195318    0.612926 B11  gam11                            :  -1.913775    1.078541 B12  gam12                            :  -1.141151    0.722682 B13  gam13                            :  -1.576021    0.844656 B14  gam14                            :  -0.565806    0.539048 B15  gam15                            :  -1.785868    0.784658 B16  gam16                            :  -0.547574    0.487481 B17  gam17                            :  -0.513223    0.522252 B18  gam18                            :  -0.942994    0.612768 B19  gam19                            :  -1.863471    0.766751       C1   eps1                             :  -3.655377    1.012951 C2   eps2                             : -28.887708  194658.590464 C3   eps3                             :  -2.574357    0.662361 C4   eps4                             :  -3.776636    1.371053 C5   eps5                             :  -2.817522    0.733333 C6   eps6                             :  -2.343885    0.644243 C7   eps7                             :  -1.969729    0.550647 C8   eps8                             :  -2.772890    0.787812 C9   eps9                             :  -1.674088    0.495354 C10  eps10                            :  -2.613472    0.741789 C11  eps11                            :  -1.806429    0.538208 C12  eps12                            :  -2.236284    0.687352 C13  eps13                            :  -1.174660    0.465809 C14  eps14                            :  -1.840413    0.622837 C15  eps15                            :  -0.991599    0.452199 C16  eps16                            :  -1.459758    0.557575 C17  eps17                            :  -1.661557    0.574279 C18  eps18                            :  -1.476907    0.515438 C19  eps19                            :  -1.139877    0.474423       D1   P[1-1]                           :   1.540445    0.212070 D2   P[2-1]                           :   1.791758    0.230284 D3   P[3-1]                           :   1.427115    0.227303 D4   P[4-1]                           :   1.930038    0.281821 D5   P[5-1]                           :   2.651756    0.354766 D6   P[6-1]                           :   1.750315    0.249536 D7   P[7-1]                           :   0.799850    0.194886 D8   P[8-1]                           :   1.195231    0.210878 D9   P[9-1]                           :   1.210567    0.208752 D10  P[10-1]                          :   1.455949    0.234064 D11  P[11-1]                          :   1.134149    0.201008 D12  P[12-1]                          :   1.207191    0.201302 D13  P[13-1]                          :   0.893969    0.209340 D14  P[14-1]                          :   0.507690    0.204052 D15  P[15-1]                          :   1.356835    0.221533 D16  P[16-1]                          :   1.555367    0.278687 D17  P[17-1]                          :   0.823684    0.210649 D18  P[18-1]                          :   1.006656    0.215069 D19  P[19-1]                          :   0.960811    0.228936 D20  P[20-1]                          :   1.324828    0.248156

 

 

*Does anyone have insight into what is happening? Am I missing something obvious?

 

Here are some other notes that might be helpful:

*of the 8,000 surveys in my data, there are 4,797 missing observations (NAs). 

*in the first year (1993), all sites were, in fact, occupied.

 

I’d be happy to send my data or attach full PRESENCE output files if it could help diagnostics.

 

Thanks!

Gavin Jones

[Jeffrey Royle]

hi Gavin,

 The first case surprises me -- i know these programs have been tested against each other in the past and they were consistent. The first model is simple enough that there should just be numerical error in the estimates I think.

 

 The 2nd example, this model probably has too many parameters per observation and it wouldn't surprise me if both solutions were local optima when they work or you have extreme sensitivity to starting values. Possibly the defaults in unmarked are not so good, not sure about PRESENCE though.

 

regards

andy

 

--
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]

A long time ago, I compared colext() vs PRESENCE vs BUGS using simulated data, and I found very similar results. However, I didn't simulate data with tons of missing values and parameters near boundaries, as you have in your data. Anyone want to do more extensive simulation tests? 

You could also try using PRESENCE's estimates as starting values in colext. This might indicate if there are real differences or if one of the programs has failed to find the global minimum. 

--

Richard Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

chandlerlab.wordpress.com

[Gavin M. Jones]

Richard and Andy--thank you for the prompt reply!

Okay, so it does appear that local minima are an issue after plugging in PRESENCE parm estimates as starting values in unmarked (seen in table below).

 

I’ve got a couple of additional questions for the group:

#1. What is going on with the AIC values? First off, the values between the two programs are in different ballparks. Additionally (and perhaps more concerning) in PRESENCE, the fully time-dependent model has all of the AIC weight—by a longshot. In unmarked, it’s the opposite (below).

 

PRESENCE

"Model"                                                      AIC deltaAIC AIC wgt Model Likelihood no.Par. -2*LogLike

"psi,gamma(year-1),eps(year-1),p(year)" 3208.34     0.00   0.9982   1.0000                      59 3090.34

"psi,gamma(.),eps(.),p(.)"                           3221.02    12.68   0.0018   0.0018                        4 3213.02

 unmarked

                                                       nPars     AIC           delta         AICwt             cumltvWt

psi(.)col(.)ext(.)p(.)                                 4      1528.04      0.00          1.0e+00         1.00

psi(.)col(year-1)ext(year-1)p(year)           59     1558.55       30.52        2.4e-07           1.00

#2. Any insight into where to go from here? I’m feeling a bit tenuous to attempt other models before I solve this basic issue. 

 

Thanks,

Gavin

 

unmarked—holding all constant

> m1 <- colext(~1, ~1, ~1, ~1, data=umf) > m1       Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,     pformula = ~1, data = umf)       Initial:  Estimate  SE      z P(>|z|)      12.3 195 0.0627    0.95   Colonization:  Estimate    SE     z  P(>|z|)     -1.66 0.246 -6.72 1.77e-11   Extinction:  Estimate    SE     z  P(>|z|)     -2.35 0.189 -12.5 1.01e-35   Detection:  Estimate     SE    z  P(>|z|)      1.21 0.0792 15.3 7.81e-53   AIC: 1528.036

unmarked—holding all constant + starting values are parm estimates from PRESENCE   > m2 <- colext(~1, ~1, ~1, ~1, data=umf, starts=c(25.064276, -0.948506, -2.07767, 1.289211)) > m2   Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,     pformula = ~1, data = umf, starts = c(25.064276, -0.948506, 2.07767, 1.289211))

Initial:  Estimate   SE      z P(>|z|)

25.1 1265 0.0198   0.984

Colonization:  Estimate    SE     z  P(>|z|)

-1.66 0.246 -6.72 1.78e-11

Extinction:  Estimate    SE     z  P(>|z|)     -2.35 0.189 -12.5 1.01e-35   Detection:  Estimate     SE    z P(>|z|)

1.21 0.0792 15.3 7.8e-53   AIC: 1528.061

 

 

unmarked—fully time dependent model   > m3 <- colext(~1, ~year-1, ~year-1, ~year, data=umf) > m3

Call: colext(psiformula = ~1, gammaformula = ~year - 1, epsilonformula = ~year -

1, pformula = ~year, data = umf)

Initial:  Estimate    SE    z  P(>|z|)

2.96 0.797 3.71 0.000205

Colonization:          Estimate      SE       z P(>|z|)

year1993 -16.6767 242.038 -0.0689  0.9451 year1994  -6.1758     NaN     NaN     NaN year1995  -8.4921     NaN     NaN     NaN year1996  -0.0866   1.545 -0.0561  0.9553 year1997  -7.6826  59.076 -0.1300  0.8965 year1998 -12.1086 208.119 -0.0582  0.9536 year1999 -10.2910 126.867 -0.0811  0.9353 year2000  -0.8808   0.853 -1.0328  0.3017 year2001  -1.3797   1.082 -1.2756  0.2021 year2002  -0.7605   0.771 -0.9861  0.3241

year2003  -1.1307   0.924 -1.2235  0.2211

year2004  -1.5021   1.137 -1.3212  0.1864 year2005  -2.0304   1.618 -1.2548  0.2096 year2006  -1.1079   0.813 -1.3634  0.1728 year2007  -1.5109   1.082 -1.3963  0.1626 year2008  -2.9904   2.703 -1.1065  0.2685 year2009  -1.8111   0.778 -2.3271  0.0200 year2010 -15.0780 220.880 -0.0683  0.9456 year2011  -2.4318   1.183 -2.0548  0.0399

Extinction:          Estimate      SE       z P(>|z|)

year1993   -11.64 113.057 -0.1030 0.91797 year1994    -3.94   2.677 -1.4725 0.14088 year1995    -9.93  56.189 -0.1768 0.85966 year1996   -17.57 241.314 -0.0728 0.94195 year1997    -3.12   1.250 -2.4939 0.01263 year1998    -1.63   0.635 -2.5610 0.01044 year1999    -1.40   0.649 -2.1651 0.03038 year2000    -2.55   1.428 -1.7824 0.07468 year2001    -1.96   0.708 -2.7655 0.00568 year2002   -13.98 207.204 -0.0675 0.94620 year2003    -1.01   0.484 -2.0810 0.03744 year2004    -1.90   0.834 -2.2769 0.02279 year2005    -2.83   2.514 -1.1240 0.26103 year2006    -2.00   0.986 -2.0248 0.04289 year2007    -2.27   0.926 -2.4558 0.01406 year2008    -1.87   0.635 -2.9444 0.00324 year2009   -12.39 136.110 -0.0911 0.92745 year2010   -15.71 230.064 -0.0683 0.94557 year2011    -1.30   0.556 -2.3456 0.01900

Detection:             Estimate    SE       z P(>|z|)

(Intercept)   0.9803 0.443  2.2152  0.0267 year1994      1.0169 0.647  1.5721  0.1159 year1995      0.1995 0.709  0.2813  0.7785 year1996      0.3429 0.618  0.5553  0.5787 year1997      1.4399 0.593  2.4288  0.0151 year1998      0.1173 0.547  0.2143  0.8303 year1999      0.0215 0.546  0.0393  0.9686 year2000      0.9807 0.667  1.4696  0.1417 year2001      0.4633 0.571  0.8108  0.4175 year2002      0.8057 0.596  1.3515  0.1765 year2003      0.1673 0.533  0.3140  0.7535 year2004      0.3709 0.526  0.7056  0.4805 year2005     -0.0181 0.583 -0.0311  0.9752 year2006     -0.8156 0.530 -1.5375  0.1242 year2007      0.4670 0.543  0.8602  0.3897 year2008     -0.1919 0.563 -0.3410  0.7331 year2009     -0.2850 0.527 -0.5412  0.5884 year2010     -0.1897 0.530 -0.3578  0.7205 year2011     -0.1741 0.585 -0.2978  0.7659 year2012      0.3507 0.564  0.6219  0.5340   AIC: 1558.554 Warning message: In sqrt(diag(vcov(obj))) : NaNs produced

unmarked—fully time dependent model + starting values are parm estimates from PRESENCE   > m4 <- colext(~1, ~year-1, ~year-1, ~year, data=umf, starts=c(StartingValues)) > m4

Call: colext(psiformula = ~1, gammaformula = ~year - 1, epsilonformula = ~year -

1, pformula = ~year, data = umf, starts = c(StartingValues))

Initial:  Estimate  SE      z P(>|z|)

24.5 419 0.0584   0.953

Colonization:          Estimate      SE         z P(>|z|)

year1993    0.198 331.081  0.000599  0.9995 year1994  -67.620 776.748 -0.087056  0.9306 year1995  -68.381 785.606 -0.087043  0.9306 year1996   -0.842   2.910 -0.289263  0.7724 year1997   -9.959 171.297 -0.058139  0.9536 year1998  -10.829 213.614 -0.050696  0.9596 year1999  -10.198 141.188 -0.072228  0.9424 year2000   -0.857   0.834 -1.027273  0.3043 year2001   -1.370   1.081 -1.267263  0.2051 year2002   -0.758   0.773 -0.980817  0.3267 year2003   -1.127   0.925 -1.218340  0.2231 year2004   -1.500   1.138 -1.318447  0.1874 year2005   -2.035   1.626 -1.251411  0.2108 year2006   -1.109   0.813 -1.363963  0.1726 year2007   -1.512   1.082 -1.397506  0.1623 year2008   -2.987   2.691 -1.110072  0.2670 year2009   -1.811   0.778 -2.327061  0.0200 year2010  -13.717 231.150 -0.059345  0.9527 year2011   -2.432   1.183 -2.054900  0.0399

Extinction:          Estimate      SE       z P(>|z|)

year1993    -3.10   1.100 -2.8183 0.00483 year1994   -28.83 447.066 -0.0645 0.94857 year1995    -9.79  35.812 -0.2733 0.78459 year1996   -13.25 136.930 -0.0968 0.92292 year1997    -3.15   1.281 -2.4615 0.01384 year1998    -1.59   0.619 -2.5727 0.01009 year1999    -1.39   0.645 -2.1609 0.03070 year2000    -2.52   1.413 -1.7836 0.07449 year2001    -1.95   0.707 -2.7582 0.00581 year2002   -11.31 110.450 -0.1024 0.91843 year2003    -1.01   0.484 -2.0768 0.03782 year2004    -1.90   0.834 -2.2770 0.02279 year2005    -2.83   2.511 -1.1255 0.26039 year2006    -2.00   0.985 -2.0259 0.04278 year2007    -2.27   0.926 -2.4561 0.01405 year2008    -1.87   0.634 -2.9439 0.00324 year2009   -12.40 152.490 -0.0813 0.93517 year2010   -13.38 209.065 -0.0640 0.94895 year2011    -1.30   0.556 -2.3455 0.01900

Detection:             Estimate    SE      z P(>|z|)

(Intercept)   0.8651 0.421  2.053 0.04012 year1994      0.9772 0.678  1.441 0.14958 year1995      0.1848 0.609  0.304 0.76142 year1996      0.4473 0.599  0.747 0.45536 year1997      1.5559 0.577  2.696 0.00703 year1998      0.2366 0.530  0.446 0.65561 year1999      0.1335 0.529  0.252 0.80073 year2000      1.1076 0.654  1.694 0.09021 year2001      0.5805 0.556  1.044 0.29630 year2002      0.9206 0.581  1.585 0.11288 year2003      0.2810 0.515  0.545 0.58547 year2004      0.4859 0.508  0.956 0.33883 year2005      0.0968 0.567  0.171 0.86456 year2006     -0.7005 0.513 -1.366 0.17208 year2007      0.5822 0.526  1.107 0.26821 year2008     -0.0769 0.546 -0.141 0.88811 year2009     -0.1696 0.509 -0.333 0.73895 year2010     -0.0744 0.513 -0.145 0.88466 year2011     -0.0589 0.569 -0.104 0.91748 year2012      0.4659 0.548  0.851 0.39483   AIC: 1560.084

To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+unsubscribe@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.

[Jeffrey Royle]

hi Gavin,

 You should never compare AIC values across software platforms because everyone leaves in or excludes certain constants from the likelihood depending on how they're building algorithms. Comparability of AIC values , therefore, should not be expected.

 

 As to whats going on with the diffeerences in estimates -- I can't imagine. I know the functions have been tested and compared in the past and we have been comfortable with it although this is making me wonder a abit

regards

andy

 

 

[Jeffrey Royle]

are these fits for a simulated data set? if not,  I think before you do anything else you should simulate some data exactly under the model you want to fit and then run it through both programs.

 

[Gavin M. Jones]

Nope, these are real data. 

I just simulated a dataset of equal dimensions and sure enough--both unmarked and PRESENCE gave nearly identical outputs for both the parameter estimates and AIC values. This was the case for both the "all constant" model and the fully time-dependent model. 

I'm going to explore this a bit further now by adding NAs into the matrix to see if the programs perform differently.