[R-sig-eco] Doubt about function colext() of package Unmarked. Wrong estimates for a metapopulation simulation.

Baldwin, Jim -FS jbaldwin at fs.fed.us
Sat Nov 24 18:34:08 CET 2012 

The problem seems to be the way you are generating the 0's and 1's and not colext in unmarked.

One can obtain a rough estimate the colonization and extinction probabilities from any two neighboring pairs of seasons (although in practice you'd want to use all of the data).  Take seasons 29 and 30:

   > table(coleta[,29],coleta[,30])

        0  1
     0 53  4
     1  4 39

An estimate for the colonization probability is 4/(4+53) = 0.07017544 and an estimate for the extinction probability is 4/(4+39) = 0.09302326.  Those are close to the colext estimates of 0.07479338 and 0.08615581.  So I would say the problem is not with colext in unmarked.

Jim


-----Original Message-----
From: r-sig-ecology-bounces at r-project.org [mailto:r-sig-ecology-bounces at r-project.org] On Behalf Of Augusto Ribas
Sent: Saturday, November 24, 2012 6:52 AM
To: r-sig-ecology at r-project.org
Subject: [R-sig-eco] Doubt about function colext() of package Unmarked. Wrong estimates for a metapopulation simulation.

Hello dear list.
I was using function colext to estimate colonization and extinction of simulated data and i'm not getting what i expected. Hope someone can enlighten me here.

The way i'm simulating is the one used in the The R Book (M Crawley), chapter 26. But its pretty straight forward. You make a scenario (matrix with 0 and 1) and loop killing populations at extinction rate then the survivors create propagulos and populate matrix.

############################################################################################
set.seed(10)
#Define Colonization (m) and extinction (e)
m<-0.20
e<-0.10


#Prepare scenario
s<-(1-e)
N<-matrix(rep(0,10000),nrow=100)
xs<-sample(1:100)
ys<-sample(1:100)
for (i in 1:100){
  N[xs[i],ys[i]]<-1
}

#Simulate 250 seasons

for (t in 1:250){
  #kill populations
  S <-matrix(runif(10000),nrow=100)
  N<-N*(S<s)
  #how many colonizaers
  im<-floor(sum( N*m))
  #place them
  placed<-matrix(sample(c(rep(1,im) ,rep (0,10000-im))),nrow=100)
  N<-N+placed
  #delete who arrived on a already colonized place
  N<-apply(N,2,function(x) ifelse(x>1,1,x)) }

#Ocuppance
sum(N)/length(N)
##################################################################################################

Then i thought i could save results from this simulation to use on package unmarked and try to figure out the colonization and extinction rates.

###############################################################################
#make a list to save results from the simulation each season
ocup<-list()

#Simulate like above and save each matrix on my ocup list for (t in 1:30){
  S <-matrix(runif(10000),nrow=100)
  N<-N*(S<s)
  im<-floor(sum( N*m))
  placed<-matrix(sample(c(rep(1,im) ,rep (0,10000-im))),nrow=100)
  N<-N+placed
  N<-apply(N,2,function(x) ifelse(x>1,1,x))
  ocup[[t]]<-N
}

#Ocupation stable
lapply(ocup,function(N) sum(N)/length(N)) ##################################################################################

Until here everything is ok. Now i try to take some samples from this population, on this 30 simulated season to use on unmarked, so i did
this:

##################################################################################

#select 100 points to sample (this is wrong right? Not all points have the same chance to be sampled)
coleta.x<-sample(1:100)
coleta.y<-sample(1:100)
data.frame(coleta.x,coleta.y)

#prepare a matrix to receive the results
coleta<-matrix(NA,ncol=30,nrow=100,dimnames=list(paste("Local",1:100),paste("Season",1:30)))
head(coleta)

for(j in 1:30) {
  for(i in 1:100) {
    coleta[i,j]<-ocup[[j]][coleta.x[i],coleta.y[i]]
  }
}

head(coleta)
##################################################################################

Then used function colext on Unmarked and:


#####################################################################################
library(unmarked)

coleta.unmarked <-unmarkedMultFrame(coleta,numPrimary=30)
coleta.prms<-colext(psiformula= ~1,gammaformula =  ~ 1,epsilonformula = ~ 1,
                    pformula = ~ 1,coleta.unmarked) coleta.prms
plogis(coleta.prms at opt$par)
################################################################################

Ang i get the results:
> coleta.prms

Call:
colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
    pformula = ~1, data = coleta.unmarked)

Initial:
 Estimate    SE     z P(>|z|)
     0.16 0.201 0.797   0.425

Colonization:
 Estimate     SE     z   P(>|z|)
    -2.52 0.0987 -25.5 2.63e-143

Extinction:
 Estimate     SE     z   P(>|z|)
    -2.36 0.0947 -24.9 3.19e-137

Detection:
 Estimate   SE     z P(>|z|)
     12.2 16.8 0.728   0.467

AIC: 1766.498


Using the logistic function to see in probability 0-1 we have:

> plogis(coleta.prms at opt$par)
[1] 0.53989597 0.07479338 0.08615581 0.99999511

I dont understand why the occupation and detectability is correct, but both colonization and extinction are not the values i used on simulation. They should be 0.20 and 0.10, not these 0.07 and 0.08.
I think maybe the way i select the sample was wrong, but anyway the answers should be near as any element of the matrix is independent.
I cant figure out why it have not worked, hope someone can easily see what i cant see here.
Does anyone know why in this example i dont get 0.20 and 0.10 as answers here(the values i used when i started the simulation up in the begin).

Well thanks for your time to read my doubt and hope you could help me figure out what i'm doing wrong here.

Have a nice weekend

Best wishes
Augusto Ribas





--
Grato
Augusto C. A. Ribas

Site Pessoal: http://recologia.wordpress.com/
Lattes: http://lattes.cnpq.br/7355685961127056

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology





This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.

More information about the R-sig-ecology mailing list