GENPRES: model with group effects?

[bglenn]

Hi -

I'm trying to conduct a power analysis for a multi-season occpuancy study where I'd like to determine the power to detect a treatment effect (10 year study, 2 groups, parameter values are the same for both groups pre-treatment, treatment occurs in year 5 - occupuancy and p change for treatment group for years 6-10-but remain constant after treatment occurs). I think what I'm trying to do is very similar to what is presented in the Bailey et al (2007) paper (Sampling design trade-offs...), where the authors estimated the power to detect a difference in occupancy between 2 groups. My problem is that I cannot figure out how to specify a model structure with group effects in GENPRES. If I generate expected values using 2 groups, how to I model that group effect (all the pre-defined models are dot or time models - do I need to create a user-specified model)? I'm sure that I'm probably missing something obvious here, but any assistance would be greatly appreciated!

Thanks,
Betsy

[darryl]

Hi Betsy,
I'm sure Jim Hines will correct me if I'm wrong, but I don't think there is a way you can currently do what you want to do within GENPRES directly. I'm pretty sure that Larissa told me that for that paper they used GENPRES to generate the data, but then took the input data, rearranged it slightly and then analyzed it manually with MARK. You could do the same by saving the input file, then manipulating the data file. Exactly what you'd have to do will depend on whether you want to use PRESENCE or MARK. Either way, it looks like Jim just appends the detection histories for each group after the previous one, so if using PRESENCE you'd just have to add to the '.pao' file a series of covariates indicating which detection history relates to each group. Or using MARK, insert an additional column and take the expected frequencies for the second group, paste them into the inserted column, then delete the original rows for the second group. Clear as mud?
Cheers
Darryl

[bglenn]

Thanks Darryl -

That does make sense (more or less:). I had considered trying to take the encounter histories into MARK, but had thought that I was just missing something obvious in GENPRES.

Betsy

[bglenn]

One more question about GENPRES - When I use a multi-season model, the encounter history generated only includes records for the primary sampling periods regardless of the number of secondary sampling periods I specify in GENPRES. If I pull these data into MARK or PRESENCE, it assumes that it is a single season dataset. Can multi season datasets with secondary sampling periods be simulated?

[jhines]

Hi Betsy,

For multi-season models, you need to specify eps,gam=0 for within-season survey intervals, and esp>0,gam>0 for between-season survey intervals. In the default scenario, GENPRES sets eps(1)=0, eps(2)=0.2, eps(3)=eps(4)=0, so it would be a 2-season dataset, with 2 surveys in the first season, and 3 surveys in the 2nd season.

As Darryl said, you can generate the data for the 2 groups in GENPRES, then create a MARK or PRESENCE input file from the GENPRES output file. Then, analyze with MARK or PRESENCE. If you're generating expected value data, this wouldn't be a bad way to go, but if you want to simulate a bunch of datasets, this would be tedious. In that case, I'd suggest writing a script in R (or whatever script language you know) to simulate the data, analyze with MARK or PRESENCE, then summarize the results. I happen to have one which was pretty close to doing what you want, so I modified it for you. Here's the R code:

Code: Select all

#
#  genpres2.R - simulate single-season data and analyze w/ PRESENCE
#                   2 groups of data are simulated: 
#                       group 1:  psi=.77,  p=.8 for all surveys
#                       group 2:  psi=.66,  p=.8 for surveys 1-2,  p=.7 for surveys 3-4
#
  #presence="h:/y/presence2/presence";
  presence="c:/progra~1/presence/presence";   #  location of PRESENCE program

  seed=-32766+floor(64542*runif(1)); TAB='\t'; #sprintf("%c",9);
  T=4; N1=500; psi1=.77; p1=c(.8,.8,.8,.8);
  T=4; N2=500; psi2=.66; p2=c(.8,.8,.7,.7);
 
  #          build design matrix input file for PRESENCE
  dm=file("genpres.dm","w");

  cat(c("0 2 3     \n"),file=dm) ;   #  MATRIX 0 - 2 rows, 3 cols
  cat(c("-,a1,a2,  \n"),file=dm,append=T)  #   matrix header line
  cat(c("psi,1,-101\n"),file=dm,append=T);  #  psi design matrix: 1=intercept, -101=1st site covariate
 
  cat(c("1",T+1,"4\n"),file=dm,append=T);  #  MATRIX 1 - T rows, 5 cols
  cat(c("-,b1,b2,b3,\n"),file=dm,append=T);            #   matrix header line
  for (i in 1:2) cat(c("p(",i,"),1,0,0\n"),file=dm,append=T,sep=""); # p part of design matrix, surveys 1-2
  for (i in 3:T) cat(c("p(",i,"),0,1,-101\n"),file=dm,append=T,sep=""); # p part of design matrix, surveys 3-4

  cat(c("2 0 0,\n"),file=dm,append=T)   #   MATRIX 2 - unused
  cat(c("3 0 0,\n"),file=dm,append=T)   #   MATRIX 3 - unused
  cat(c("4 0 0,\n"),file=dm,append=T)   #   MATRIX 4 - unused
  cat(c("5 0 0,\n"),file=dm,append=T)   #   MATRIX 5 - unused
  close(dm)
   
  #    generate data for group 1
  present1=(runif(N1)<psi1); det1=rep(0,N1*T); dim(det1)=c(N1,T); det2=det1;
  for (i in 1:T) det1[,i]=(runif(N1)<p1[i]);
  present2=(runif(N2)<psi2);     #    generate data for group 2
  for (i in 1:T) det2[,i]=(runif(N2)<p2[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+N2,"\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,present2*det2),file=pao,append=T,sep='\t',row.names=F,col.names=F); # output detection data
  #     write number of site covariates, followed by the name of the 1st (and only) site covariate...
  cat("1\ngrp2\n",file=pao,append=T,sep="");
  #     write site covariate vector...
  write.table(rbind(rep(0,N1),rep(1,N2)),file=pao,append=T,sep='\n',row.names=F,col.names=F); # output detection data
  #     write number of survey covariates (zero) and close input file...
  cat("0\n",file=pao,append=T,sep="");
  close(pao)
  #                             run PRESENCE... args: i=inputfilename  model=100 for single-season model
  cmd=sprintf("%s i=genpres.pao model=100 name=psi_p QUIET",presence)
  system(cmd)
  #     scan PRESENCE output file and save parameter estimates...
  v=readLines("presence.out"); l=length(v); est="";
  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[4],c[5]);       #  so we can save the 4th and 5th fields (the param estimate and std.err)
    }

  cat(c(">>",est,"\n"),sep="\n");  #   print estimates


[bglenn]

Thanks Jim-

That helps a lot (especially the code). GENPRES doesn't seem to run if > 15 survey periods are specified, so I think that for what I'm doing a script is the way to go.

Betsy

[jhines]

Betsy,

If you have a lot of surveys, GENPRES will take a long time to run, unless you check the option 'simulate by individual'. When this option is checked, the process draws random numbers for each site to determine occupancy, then draws random numbers to determine if the species is detected on each survey. So, it only has to loop through the individual sites, and surveys. When it is not checked, GENPRES will try to keep track of every single possible detection history. If you have 15 surveys, there are 2^15 (=32768) possible histories.

[bglenn]

Thanks again - That makes sense (now that I know what 'simulate by individual' does).