simulations of gall damage

galls.R

## a simulation of detecting effect size on gall parasitisim rates
## and occupancy.
## basically ask: how many larvae are in each gall (poisson)
## how many samples do you need to see the difference?
##

#first, a homemade function for the Zero-Inflated Negative Binomial distribution
rzinbinom <- function(n,mu,size,zprob){
  ifelse(runif(n)<zprob,0,rnbinom(n,mu=mu,size=size))
}

n <- 10  #the number of genotypes we're studying
salix <- expand.grid(replicate=1:12,genotype=letters[1:n]) # the beginning of the dataframe: replicates of each genotype

mean.galls <- round(runif(n=n,min=50,max=500)) #the mean number of galls per genotype
gall.fail <- runif(n) #the probability of an aborted gall for each genotype
mean.per.gall <- round(runif(n=n,min=1,max=10))#the mean density of insects per gall for each genotype

salix <- within(salix,{
  tree <- paste(genotype,replicate,sep="")
  gall.density <- sapply(genotype,function(x) rpois(1,mean.galls[as.numeric(x)])) #generate gall density per rep
  prob.gall.fail <- abs(sapply(genotype,function(x) rnorm(1,gall.fail[as.numeric(x)],sd=0.1)))
  prob.gall.fail[which(prob.gall.fail>1)] <- 1-(prob.gall.fail[which(prob.gall.fail>1)]-1)
  n.per.gall.true <- sapply(genotype,function(x) rpois(1,mean.per.gall[as.numeric(x)]))
})

#with(salix,plot(prob.gall.fail~genotype))
with(salix,plot(n.per.gall.true~genotype))
with(salix,plot(gall.density~genotype))

all.galls.rbinom <- with(salix,mapply(rzinbinom,
                                      n=gall.density,
                                      mu=n.per.gall.true,
                                      size=3,
                                      zprob=prob.gall.fail
                                      )
                         )
                                      

sampled.galls <- lapply(all.galls.rbinom,function(x) sample(x,size=10,replace=FALSE))
galls <- data.frame(do.call(rbind,sampled.galls))
dimnames(galls)[[2]] <- paste("g.",1:n,sep="")
salix$occ <- rowSums(galls>0)/n
salix$n.per.gall <- apply(galls,1,function(x) mean(x[x>0])
                          )

par(mfrow=c(1,2))
#with(salix,plot(occ~genotype))
with(salix,plot(occ~prob.gall.fail))
#with(salix,plot(n.per.gall~genotype))
with(salix,plot(n.per.gall~n.per.gall.true))
abline(a=0,b=1)

tapply(salix$n.per.gall,salix$genotype,mean)