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

[Augusto Ribas]

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