creating the unmarkedMultFrame

[Joaquín Aldabe]

Dear all, this is Joaquín Aldabe, ecologist and ornithologist from Uruguay. I´m new in this forum, and I look forward to learning and, if possible, help others. I´m trying to creat a unmarkedMultFrame for running a dynamic occupancy model. 

I´ve been reading the R help and the article by Kery and Chandler (2016) but still have doubts. I´m confused about the difference between siteCovs, obsCovs and yearlySiteCovs. 

My data set: two years where I visited 81 sites; sites were visited between 1 and 4 times (for probability of detection estimation). I have three covariates that vary among sites and between the two years, and I wanto to assess if these covariates govern colonization and extinction rates. Covariables are grass height, seed abundandance (food for Marsh Seedeater) and number of perches. I also have a variable that may affect detection probability (wind intensity, that usually affect birds activity and hence probability of detection).

In year 1, I did not do repeated visits, but will assume probability of detection si similar to year 2 (that I do have repeated surveys). 

What covariate should have siteCovs, obsCovs and yearlySiteCovs respectively.

Thanks for your time in advanced.

Best regards,

Joaquín. 

[Kery Marc]

Dear Joaquin,

 

You will have to organize detection data in a “balanced” way, i.e., fill up the missing occasions in year 1 with NAs. Hence, the top of your detection data, as you must supply it for unmarked, will look like this:

 

1-NA-NA-NA-0-1-0

0-NA-NA-NA-0-0-0

.....

Then, by telling unmarked that you have two years and 4 occasions each, things will be filled in in the right way. You should then NOT fit models where p is made year-dependent, because you must make p identical across years in order to inform p in year 1 by the information coming from year 2.

 

As for the types of covariates, they can be somewhat confusing perhaps when you only have two occasions. In this case, one of the three types doesn’t make sense: you cannot fit yearly site-covariates in a study with only 2 primarly occasions, since yearly site-covariates are used to explain variation in the dynamic rate params (colonisation, extinction). But with only two primary occasions, you only have a single colonisation and a single extinction parameter, hence, you can only try to explain spatial, but not spatio-temporal variation in them. That’s done by fitting site covariates.

 

Best regards  --- Marc

--
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/d/optout.

[Joaquín Aldabe]

Ok thankyou very much Marc! Now it worked the unmarkedMultFrame. So you mean that I should model occupancy with "occu" function? and not estimate colonization and extinction parameter with colex function?

When you say 

" But with only two primary occasions, you only have a single colonisation and a single extinction parameter, hence, you can only try to explain spatial, but not spatio-temporal variation in them", you mean in colonisation and extinction parameters? I tried this doing 

fm2<-colext(psiformula=~y1.grass+y1.seed+y1.perch,

            gammaformula=~y2.grass+y2.seed+y2.perch,

            epsilonformula = ~y2.grass+y2.seed+y2.perch,

            pformula=~1,data=umf1)

But it did not work. y1.grass is grass height in the first year, y2. grass is the grass height in the second year, and so on with seed and perch variables. 

Cheers,

Joaqu'in. 

[Kery Marc]

Dear Joaquin,

glad you get stuff into an unmarkedMultFrame.

But no, I was not suggesting to model your problem with occ. I think that it can be interesting to partition the change in occupancy between period 1 and 2 into the two components: colonisation and extinction. These are different processes and may be affected by different covariates.

But if you are not interested in these processes, you can simply stack the data, treat year as a factor and use occu to fit a static model. That's an alternative.

If you want to model colonisation/extinction dynamics then you must note that you cannot fit a yearly site-covariate into the colo or the ext models, only site covariates work. With two primary periods, you only have a single parameter each for colo and extinction and hence all you can model there is spatial heterogeneity. This is taken account for by a site covariate.

Re. those covariates that you have measured at both primary occasions: in the context of a model with colo and extinction I would think about what you know about these birds mechanistically. Is it the environmental conditions during the first or during the second primary occasion that seems more likely to affect colo or extinction ? Then take that. If you can't decide, then you could take the average and use that as a site covariate. Or, you specify these covariates for each year as a separate site covariate and fit each of them. But check first their correlation: if they're highly correlated (e.g., r > 0.6 or 0.7) then there is no use fitting them both.

Best regards  --- Marc

---

From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Joaquín Aldabe [joaquin...@gmail.com]
Sent: 22 August 2017 18:59
To: unmarked
Subject: [unmarked] Re: creating the unmarkedMultFrame

[Joaquín Aldabe]

Hi Marc, I want to model colonisation/extinction dynamics. I think grass height, as well as seed abundance and number of perches are determining first year occupancy. As for colonization, I think it is positively affected by the increase of the number of perches per site in the second year. Extinction may be influenced by a decrease in seed abundance in the second year. So, I typed this script, which does not work (it doesn't create the object):

fm1 <- colext(psiformula = ~y1.seed+y1.perch, # First-year occupancy

             gammaformula = ~y2.perch , # Colonization

             epsilonformula = ~ y2.seed, # Extinction

             pformula = ~ 1, # Detection

             data = umf1)

Here is str of the unmarkedMultFrame. Is it correct?

> str(umf1)

Formal class 'unmarkedMultFrame' [package "unmarked"] with 7 slots

  ..@ numPrimary    : num 2

  ..@ yearlySiteCovs: NULL

  ..@ y             : int [1:81, 1:8] 1 1 1 1 1 1 1 1 1 1 ...

  .. ..- attr(*, "dimnames")=List of 2

  .. .. ..$ : NULL

  .. .. ..$ : chr [1:8] "y1.v1" "y1.v2" "y1.v3" "y1.v4" ...

  ..@ obsCovs       : NULL

  ..@ siteCovs      :'data.frame': 81 obs. of  6 variables:

  .. ..$ y1.grass: int [1:81] 74 70 74 73 121 200 138 90 86 45 ...

  .. ..$ y2.grass: num [1:81] 82 NA 67 102 64 160 100 130 112 180 ...

  .. ..$ y1.seed : int [1:81] 77 4 59 175 50 72 85 109 77 39 ...

  .. ..$ y2.seed : int [1:81] 30 NA 65 43 14 60 52 45 40 0 ...

  .. ..$ y1.perch: int [1:81] 28 31 123 116 20 25 34 33 35 18 ...

  .. ..$ y2.perch: int [1:81] 3 NA 54 60 5 1 4 3 19 12 ...

  ..@ mapInfo       : NULL

  ..@ obsToY        : num [1:8, 1:8] 1 0 0 0 0 0 0 0 0 1 ...

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/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/SnGhBynqJPU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Joaquín Aldabe

Grupo Biodiversidad, Ambiente y Sociedad
Centro Universitario de la Región Este, Universidad de la República

Ruta 15 (y Ruta 9), Km 28.500, Departamento de Rocha  

                                                           

Departamento de Conservación
Aves Uruguay
BirdLife International
Canelones 1164, Montevideo

https://sites.google.com/site/joaquin.aldabe

[Kery Marc]

Dear Joaquin,

what is the error message you get when trying to fit that model ? And can you show us summary(umf1) ?

Best regards  --- Marc

[Joaquín Aldabe]

This is the error message:

Error in optim(starts, nll, method = method, hessian = getHessian, ...) :
  initial value in 'vmmin' is not finite
In addition: Warning message:
Some observations have been discarded because correspoding covariates were missing.

This is the summary of the umf1 (it is also attached in Tinn R format as it seems that in the body text the columns move):

unmarkedFrame Object

81 sites
Maximum number of observations per site: 8
Mean number of observations per site: 3.59
Number of primary survey periods: 2
Number of secondary survey periods: 4
Sites with at least one detection: 34

Tabulation of y observations:
   0    1 <NA>
 233   58  357

Site-level covariates:
    y1.grass         y2.grass         y1.seed         y2.seed          y1.perch         y2.perch    
 Min.   : 15.00   Min.   :  0.00   Min.   :  0.0   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
 1st Qu.: 52.00   1st Qu.: 52.65   1st Qu.:  2.0   1st Qu.:  0.00   1st Qu.:  1.00   1st Qu.:  3.00  
 Median : 75.00   Median : 80.00   Median : 18.0   Median : 10.00   Median : 12.00   Median :  9.50  
 Mean   : 91.62   Mean   : 95.40   Mean   : 32.6   Mean   : 20.50   Mean   : 25.64   Mean   : 34.59  
 3rd Qu.:108.00   3rd Qu.:120.00   3rd Qu.: 50.0   3rd Qu.: 38.25   3rd Qu.: 32.00   3rd Qu.: 33.75  
 Max.   :280.00   Max.   :240.00   Max.   :177.0   Max.   :100.00   Max.   :200.00   Max.   :556.00  
                  NA's   :2                        NA's   :3                         NA's   :3    

[Kery Marc]

Hi Joaquin,

OK, so nothing totally wrong, except that the function optimiser can't find a solution. Try this: start with the simplest model with only intercepts and then progress in small steps, by adding one covariate at a time and using the solutions of the previous model as starting values for the current model fit.

Best 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
______________________________________________________________

***  Hierarchical modeling books   ***
(4) Kéry & Royle (2017): AHM Volume 2,  Dynamic and advanced models, in prep.
(3) Kéry & Royle (2016): Applied hierarchical modeling in ecology, Academic Press, Volume 1, Prelude and Static Models, see www.mbr-pwrc.usgs.gov/pubanalysis/keryroylebook/
(2) Kéry & Schaub (2012): Bayesian Population Analysis using WinBUGS, Academic Press; see www.vogelwarte.ch/bpa
(1) Kéry (2010): Introduction to WinBUGS for Ecologists, Academic Press; see www.mbr-pwrc.usgs.gov/pubanalysis/kerybook 

***   Hierarchical modeling workshops ***
www.phidot.org/forum/viewforum.php?f=8

***   Hierarchical modeling Google Group mailing lists   ***
(1) unmarked: for questions specific to the R package unmarked
(2) hmecology: for material covered in Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2015); see groups.google.com/forum/?hl=en#!forum/hmecology

---

From: Joaquín Aldabe [joaquin...@gmail.com]
Sent: 23 August 2017 15:39
To: Kery Marc
Cc: unma...@googlegroups.com

[Joaquín Aldabe]

It worked when I modeled psiformula with covariates (scaled and centered):

model.1 <- colext(psiformula = ~y1.seed+y1.perch, # First-year occupancy

             gammaformula = ~1 , # Colonization

             epsilonformula = ~ 1, # Extinction

             pformula = ~ 1, # Detection

             data = umf1)

But when I add even just one variable to gammaformula or epsilonformula it shows me the same error. For example:

model.1 <- colext(psiformula = ~1, # First-year occupancy

              gammaformula = ~1 , # Colonization

              epsilonformula = ~ y2.seed, # Extinction

              pformula = ~ 1, # Detection

              data = umf1)

Any idea?

[Kery Marc]

Perhaps you don't have very many observed state change events ? Can you show us summary(model.1) ?

---

From: Joaquín Aldabe [joaquin...@gmail.com]
Sent: 23 August 2017 16:17

[Joaquín Aldabe]

sure. Here it is:

summary(model.1)

Call:

colext(psiformula = ~y1.seed + y1.perch, gammaformula = ~1, epsilonformula = ~1, 

    pformula = ~1, data = umf1)

Initial (logit-scale):

            Estimate    SE     z P(>|z|)

(Intercept)    1.167 1.030 1.133  0.2571

y1.seed        4.135 1.789 2.311  0.0208

y1.perch       0.535 0.741 0.723  0.4699

Colonization (logit-scale):

 Estimate    SE     z P(>|z|)

    -1.28 0.504 -2.55  0.0108

Extinction (logit-scale):

 Estimate    SE     z P(>|z|)

    0.185 0.471 0.392   0.695

Detection (logit-scale):

 Estimate    SE      z P(>|z|)

  -0.0337 0.237 -0.142   0.887

AIC: 257.0085 

Number of sites: 81

optim convergence code: 0

optim iterations: 38 

Bootstrap iterations: 0 

[Kery Marc]

don't see anything conspicuous. You should try to provide starting values then.

This goes something like this:
(1) get the MLEs for the simpler model, e.g. like this: coef(model.1)
(2) assume that you want to fit a model where you have one additional covariate in the model for colonisation. Then, you have to squeeze into the right position of the vector of coefficients a zero and provide that vector for starting values, using the starts argument in the colext fitting function.

Best regards  --- Marc

---

From: Joaquín Aldabe [joaquin...@gmail.com]
Sent: 23 August 2017 16:29

[Joaquín Aldabe]

ok. This is what I did. 

ran coef function for getting MLE of the coefficients of the model with covariates associated just to psi:

coef(model.1)

     psi(Int)  psi(y1.seed) psi(y1.perch)      col(Int)      ext(Int) 

   1.16699470    4.13455875    0.53512199   -1.28441477    0.18498429 

       p(Int) 

  -0.03367617

Then, I made a vector with each MLE value, and added a cero to the right of the vector. Like this: 

coef.0<-c(1.16699470,4.13455875,0.53512199,-1.28441477,0.18498429,-0.03367617,0)

Then I runned 

model.2 <- colext(psiformula = ~y1.seed+y1.perch, # First-year occupancy

                  gammaformula = ~y2.perch , # Colonization

                  epsilonformula = ~ 1, # Extinction

                  pformula = ~ 1, # Detection

                  data = umf1, starts = coef.0)

But had the same error..

Am I correctly following your instructions?

Thanks so much, joaquin. 

[Kery Marc]

Dear Joaquin,

no, you don't fill in the 0 at the right place in the vector of inits. You must provide the inits in exactly the same order as the coefficients are ordered in the model. The coefficient for y2.perch in the submodel for colonization comes after the colonization intercept. Therefore, you must provide the starting values in this order (note the filled in 0 in red).

starts <-c(1.16699470,4.13455875,0.53512199,-1.28441477, 0, 0.18498429,-0.03367617)

Then, you do

model.2 <- colext(psiformula = ~y1.seed+y1.perch, # First-year occupancy

                  gammaformula = ~y2.perch , # Colonization

                  epsilonformula = ~ 1, # Extinction

                  pformula = ~ 1, # Detection

                  data = umf1, starts = starts)

Best regards  --- Marc

---

From: Joaquín Aldabe [joaquin...@gmail.com]
Sent: 23 August 2017 17:28