Hi!
I have a question:
Would you like to programming the Presence in a R packages, like Density (secr)?
Thanks
#
# genpres3.R - simulate single-season data and analyze w/ PRESENCE
# psi=.77, p=.8 for all surveys
#
presence="c:/progra~1/presence/presence"; # location of PRESENCE program
seed=-32766+floor(64542*runif(1));
T=8; N1=500; psi1=.77; p1=rep(.6,T); nsims=2;
# build design matrix input file for PRESENCE
dm=file("genpres.dm","w");
cat(c("0 4 4 \n"),file=dm) ; # MATRIX 0 - 3 rows, 3 cols (1 extra row,col for labels)
cat(c("-,a1,a2,a3 \n"),file=dm,append=TRUE) # matrix header line
cat(c( "psi,1,0,0\n"),file=dm,append=TRUE); # psi design matrix:
cat(c("ssPsi1,0,1,0\n"),file=dm,append=TRUE); # ssPsi1 (alias: theta0)
cat(c("ssPsi2,0,0,1\n"),file=dm,append=TRUE); # ssPsi2 (alias: theta1)
modl=1; # set modl=1 for constant detection (p(.)), modl=2 for p(t)
if (modl==1) {
cat(c("1",T+1,"2\n"),file=dm,append=TRUE); # MATRIX 1 - T rows, 1 cols
cat(c("-,b1\n"),file=dm,append=TRUE); # matrix header line
# write design matrix for p ... model: intercept only ... T rows, 1 col
write.table(cbind(paste("p",1:T,sep=""),rep(1,T)),row.names=F,col.names=F,quote=F,sep=",",file=dm,append=TRUE)
} else {
cat(c("1",T+1,"T+1\n"),file=dm,append=TRUE); # MATRIX 1 - T rows, T cols
cat(c("-",paste("b",1:T,sep=""),"\n"),file=dm,append=TRUE); # matrix header line
# write design matrix for p ... model: p(t) ... T rows, T cols
write.table(cbind(paste("p",1:T,sep=""),diag(T)),row.names=F,col.names=F,quote=F,sep=",",file=dm,append=TRUE)
}
cat(c("2 0 0,\n"),file=dm,append=TRUE) # MATRIX 2 - unused
cat(c("3 0 0,\n"),file=dm,append=TRUE) # MATRIX 3 - unused
cat(c("4 0 0,\n"),file=dm,append=TRUE) # MATRIX 4 - unused
cat(c("5 0 0,\n"),file=dm,append=TRUE) # MATRIX 5 - unused
close(dm)
if (file.exists("presence.out")) file.remove("presence.out");
if (file.exists("genpres_sims.csv")) file.remove("genpres_sims.csv");
for (isim in 1:nsims) { cat(c("simulation ",isim,"\n"));
# generate data for group 1
present1=(runif(N1)<psi1); det1=rep(0,N1*T); dim(det1)=c(N1,T);
for (i in 1:T) det1[,i]=(runif(N1)<p1[i]);
# create PRESENCE input file: genpres.pao
# 1st line in input file contains N and T, separated by TAB char,
# where N=total number of sites, T=number of surveys
pao=file("genpres.pao","w")
cat(c(N1,"\t",T,"\n",sep=""),file=pao,append=F,sep=""); # output nsites,nsurveys
# write detection histories to input file...
# the two groups of site-histories are concatenated, and the detections are
# multiplied by the random presence variable
write.table(rbind(present1*det1),file=pao,append=TRUE,sep='\t',row.names=F,col.names=F); # output detection data
# write number of site covariates (zero)
cat("0\n",file=pao,append=TRUE,sep="");
# write number of survey covariates (zero)
cat("0\n",file=pao,append=TRUE,sep="");
# write number of surveys per season and close input file...
cat(c(T,"\n"),file=pao,append=TRUE,sep="");
close(pao)
# run PRESENCE... args: i=inputfilename
# model=1000 for single-season model
# model=7000 for single-season, spatial-dep. model
# design matrix file assumed to be same as input file with extention 'dm' replacing 'pao'
cmd=sprintf("%s i=genpres.pao model=7000 name=psi_th_p QUIET",presence)
system(cmd,ignore.stdout=TRUE,ignore.stderr=TRUE)
# scan PRESENCE output file and save parameter estimates...
v=readLines("presence.out"); l=length(v); est=NULL;
ii=grep("survey_",v); l=length(ii);
for (j in 1:l) {
a=v[ii[j]]; # a=jth line of output file
b=gsub(" +","@",a); # replace multiple blanks with a single @ char
c=unlist(strsplit(b,"@")); # split output line into fields
est=c(est,c[5],c[6]); # so we can save the 5th and 6th fields (the param estimate and std.err)
}
dim(est)=c(2,length(est)/2);
#print(t(est),quote=F,col.labels=c("estimate","std.err")); # print estimates
cat(c(est,"\n"),sep=",",file="genpres_sims.csv",append=TRUE); # save estimates to csv file
}
Hi,
First of all, thank you for the code and the quick answer!
I understand how it workes, but can't get an .out file, thus the presence can not estimate the parameters!
I tried on windows 7, XP with Presence v3.1 <110411.1002>
cmd=sprintf("%s i=genpres.pao model=7000 name=psi_th_p QUIET",presence)
system(cmd,ignore.stdout=TRUE,ignore.stderr=TRUE)
system(paste('"c:/Program Files (x86)/Presence/presence.exe"','i=genpres.pao model=7000 name=buco'),ignore.stdout=TRUE,ignore.stderr=TRUE)
Return to software problems/news
Users browsing this forum: No registered users and 11 guests