R code for single season single species models?

[R.S.]

Hi everyone; are R codes for the single season single species models available? I was urged by my supervisor to move away from the comfi Presence interface ...

[jlaake]

I've seen a message that the singe season occumapncy model iis in the package called unmarked on CRAN for R. I've not used it. The single season occupancy model can also be run in MARK through RMark.

What is the reason for not using Presence in this case? Are you doing simulation or are you trying to script the analysis? Or does your adviser just want you to learn something different?

[R.S.]

Thanks! The reason for not using Presence is that I have to learn R eventually, and better sooner than later - for some planned analyses - and that I might as well do it with my own data and these analyses... It's a little scary, though.

[cnagy]

these helped me to get started

http://www.stat.washington.edu/cggreen/rprimer/

and this site for specific analyses

http://www.ats.ucla.edu/stat/R/dae/default.htm

[jlaake]

Good for your adviser! Offlist I'll send you a 50+ page mini-tutorial of sorts that I put together that covers some of the basic aspects of R. There is lots of free material on the web to help you learn R, Look under Documentation on the R Home page. I'd suggest one or more books on R as well. One I suggest is Data Manipulation in R by P. Spector. There are numerous Use R books available via Springer. If you can only afford one book then you may want to consider the R Book by Crawley published by Wiley. It is quite comprehensive but don't use attach like he does in the book. If you have never had any programming then it may be a bit daunting at first but in the long run I think it is good advice.

[R.S.]

Thanks a lot, folks. Truly appreciate any advise or hint!

[jhines]

Here's some R code for the simplest single-season model using sample data included in PRESENCE:

Code: Select all

#   single-season occupancy model using Blue Ridge Salamander data
a=c(0,0,0,1,1,
0,1,0,0,0,
0,1,0,0,0,
1,1,1,1,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
1,0,0,0,0,
0,0,1,1,1,
0,0,1,1,1,
1,0,0,1,1,
1,0,1,1,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,1,0,
0,0,0,1,0,
0,0,0,0,1,
0,0,0,0,1);
n=length(a)/5; dim(a)=c(5,n); b=t(a); eps=.0000000000000001;

lik <- function(th) {      #   likelihood function - th=parameter vector
  lik=0; psi=1/(1+exp(-th[1])); p=c(0,0);
  p[2]=1/(1+exp(-th[2])); p[1]=1-p[2];
  for (i in 1:n) {           #  loop through n sites
    prd=psi*prod(p[1+b[i,]]);   # det. hist prob. = product(psi*pq(j))
    if (sum(b[i,])<1) prd=prd+1-psi;  #  where pq(j)=p if detected in survey j
    lik=lik+log(prd)                  #        pq(j)=1-p if not detected in j
  }
  return(-lik)
}
out <- nlm(lik,c(0.1,0.2),hessian=T,print.level=0) # find max likelihood value
beta=out$estimate; vc=solve(out$hessian);           #  = min neg. likelihood value
print(out); cat("SE(theta)="); print(sqrt(diag(vc)))
psi=1/(1+exp(-beta[1])); p=1/(1+exp(-beta[2]));
d=exp(-beta[1])*psi*psi; sepsi=sqrt(d*vc[1,1]*d); # get se(psi) via delta method
d=exp(-beta[2])*p*p; sep=sqrt(d*vc[2,2]*d);      # 
cat(c("psi=",psi," se=",sepsi,"\n"));
cat(c("p=",p," se=",sep,"\n"));



[bacollier]

Looks like everyone is getting in on the action, so here is an example in RMark.

R.S. Here is a RMark example using the data from RDOccupancy in RMark without the occasion specific information (e.g., just the ch, samplearea, and cover columns were used). Note you will have to use RMark (see the Appendix for specifics) and change your working directories and such but this should get you started.


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Example application of RMARK for single season occupancy modeling-ROccupSingleSeason_example #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#load RMark
library(RMark)

# set working directory for analysis dataset
setwd("C:/Documents and Settings/bacollier.ECOSYSTEM/Desktop/TX A&M University Projects/RmarkRDExample/SingleSeasonOcc/") # will need to change to your location

# read in example encounter history which is a simple tab-delimited text file that looks
# like a MARK .inp file but (1) does not contain any comments or semicolons and (2) includes a 1st row
# that contains the column labels, and (3) does not contain a column designating unique encounter histories

ROccup<-read.table("ROccupSingleSeason_example.txt", colClasses=c("character", rep("numeric", 2)), header=TRUE)
ROccup

run.ROccup=function()
{
psi.dot=list(formula=~1)
psi.cover=list(formula=~cover)
psi.samplearea=list(formula=~samplearea)

p.dot=list(formula=~1)
p.time=list(formula=~time)

ROccup.process=process.data(ROccup, model="Occupancy")
ROccup.ddl=make.design.data(ROccup.process)

#Candidate Model List
model.1<-mark(ROccup.process, ROccup.ddl, model.parameters=list(p=p.dot, Psi=psi.dot), invisible=FALSE)
model.2<-mark(ROccup.process, ROccup.ddl, model.parameters=list(p=p.time, Psi=psi.cover), invisible=FALSE)
return(collect.models())
}


ROccup.example<-run.ROccup() #This runs the above models

#Allows you to open the standard MARK .txt output file

ROccup.example$model.1

[darryl]

Yet another option, is that you can call PRESENCE from within R (similar to how RMark interfaces with Mark), you just need to set up design matrices etc within R. I think Jim Hines has posted something on that before. It all depends on how much of the wheel you (or your supervisor) want to reinvent.

Cheers
Darryl

[jlaake]

At the risk of turning the phidot forum into an R forum, below I'm posting a revision to Jim's single season occupancy model. But first let me say, I'm glad to see Jim jumping into R and I learned something from his posting. When I learned programming it was in FORTRAN probably like Jim did. In switching from something like FORTRAN to R it is often hard to adjust to vectorized operations in R. A few years ago I would have likely posted a similar solution that Jim did. As I've learned more about R I've learned to adopt vectorized thinking. This has a couple of advantages. First it tends to speed up the code although these days loops in R are appreciably faster than they used to be. But more importantly, it tends to make solutions more generalizable. I took Jim's lik function, changed the computation to use a vectorized form with a matrix of p's (sites by time) and a vector of Psi's (by site) and the code can now be used to solve any single season occupancy model with formula rather than just the specific p(.)Psi(.) model that Jim's would fit. I didn't have to create a DM by hand and I didn't use any concept of a PIM. That's because I created a dataframe containing a record for each site-occasion pairing. If MARK was written that way there would be no need for PIMS. To be honest, I've not used PRESENCE so I don't know if it has PIMS. PIMS and creating DMs by hand needs to go away. It was a functional concept at one point in time but we should be moving away from it. If some of you caught Evan's presentation of our talk at EURING, he mentioned this at the end. Anyhow, I'll stop my blathering. Below is the code with the same data and some added "habitat" data to make it interesting. To make this generally useful, you would have to write additional code to add variable names to estimates and have a summary etc. No need to do that because the R package unmarked should do that, although I've not used it to be sure.

a=c(0,0,0,1,1,
0,1,0,0,0,
0,1,0,0,0,
1,1,1,1,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
0,0,1,0,0,
1,0,0,0,0,
0,0,1,1,1,
0,0,1,1,1,
1,0,0,1,1,
1,0,1,1,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,0,0,
0,0,0,1,0,
0,0,0,1,0,
0,0,0,0,1,
0,0,0,0,1);
k=5
n=length(a)/k; dim(a)=c(k,n); b=t(a); eps=.0000000000000001;
occdata=data.frame(Site= rep(1:n,each=k),habitat=rep(c(rep("a",20),rep("b",19)),each=k),time=factor(rep(1:k,n)))

seen=as.numeric(rowSums(b)>0)
lik <- function(th) { # likelihood function - th=parameter vector
psi=plogis(DM.psi%*%th[1:ncol(DM.psi)])
p=matrix(plogis(DM.p%*%th[(ncol(DM.psi)+1):length(th)]),nrow=n,ncol=k,byrow=TRUE)
neg.loglik.occ=-sum(seen*(b*log(p)+(1-b)*log(1-p)))-sum(seen*log(psi))
neg.loglik.notocc=-sum((1-seen)*log(apply((1-p)^(1-b),1,prod)*psi+1-psi))
return(neg.loglik.occ+neg.loglik.notocc)
}
DM.p=model.matrix(~1,occdata)
DM.psi=model.matrix(~1,occdata[occdata$time==1,])
out.null <- nlm(lik,c(0.1,0.2),hessian=T,print.level=0)

DM.p=model.matrix(~time,occdata)
DM.psi=model.matrix(~habitat,occdata[occdata$time==1,])
out.time.hab <- nlm(lik,c(0.1,0,0.2,0,0,0,0),hessian=T,print.level=0)