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-boun...@r-project.org
[mailto:r-sig-ecology-boun...@r-project.org] On Behalf Of Augusto Ribas
Sent: Saturday, November 24, 2012 6:52 AM
To: r-sig-ecology@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@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@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@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.
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
