ип# Program to simulate krill sex ratio # and carry out sensitivity to assumptions # # S.G. Candy # Vers 2.0 (15 Jan 2007) # # method: # # 1. Intialize input parameters : # * Number of individuals of each sex and age (before mortality imposed)in a simulation # * VB parameters for each of Male and Female # * number of age classes # * mortality rate parameter, M, for females and ratio of Male to Female mortality rate # * standard deviation of lognormal errors for sex-specific values of M # * standard deviations for each length about the mean for each age class and sex # Note: Load functions "simulate.krill.call" and "simulate.krill.sexr" at the end of this script # # 2. Call function to carry out steps 3 to 5 # # 3. assign to 10,000 # individuals in each age class 1 to A and each sex a random normally # distributed length about the VB predicted mean using a standard deviation # inputed for each age class # # 4. Randomly remove individuals from each sex and age class that have 'died' # using an instantaneous mortality rate of M that can vary between sexes # and can be accelerated in the males using a value of lambda > 0 # # 5. Accumulate numbers for each sex across ages in each 1 mm length bin from # 20 to 70 mm and calculate the sex ratio for each bin and graph # # 6. Repeat steps 3 to 5 to incorporate uncertainty in parameter values stochastically by drawing parameter # values from distributions of parameterA values (only M's at this time) about mean values used in (4) # # 7. Repeat steps 1 to 6 with different fixed parameter values # # ##################################### # Step 1 # ##################################### library(lattice) No.ages.m <- 7 No.ages.f <- 7 # mortality M.mr is the ratio of male to female M parameter #M.f <- 0.66 M.f <- 0.8 M.mr <- 1.0 #M.mr <- 0.66 M.m <- M.mr*M.f # Input lambda parameter value for exponential departure male mortality # (i.e.lambda>0, lambda=0 corresponds to no acceleration after age age0) age0 <- 3.0 lambda <- 1.5 # input Standard deviation for lognormal errors in M.f if stochastic variation required # Mort.m.sd > 0 allows stochastic variation in the relationship between M.f and M.m Mort.sd <- 0.2 Mort.m.sd <- 0.2 #Mort.sd <- 0.0 #Mort.m.sd <- 0.0 # No.ind is number of individuals in a simulation # No.run is the number of runs of the simulation in which Mp.f and Mp.m # can be varied stochastically No.ind <- 10000 No.run <- 1000 # ages are assumed to be the sequence 1 to No.ages #VB.m <- c(66.84963,0.43747,-0.22435) #VB.f <- c(54.63266,0.601248,-0.19512) VB.m <- c(60.1334,0.4981,-0.3511) VB.f <- c(60.1334,0.4308,-0.3511) #VB.f <- VB.m # age-class standard deviations (there should be one for each age class) SD.m <- c(2.0,2.5,3.0,3.5,4.0,4.5,5.0) SD.f <- c(2.0,2.5,3.0,3.5,4.0,4.5,5.0) # means by age classes Avlen.m <- VB.m [1] * (1 - exp(-VB.m[2]*(c(1:No.ages.m)-VB.m[3]))) Avlen.f <- VB.f [1] * (1 - exp(-VB.f[2]*(c(1:No.ages.f)-VB.f[3]))) # set up vector of 3 quantiles to summarise data # eg. median, lower 5%, upper 95% quant <- c(0.5,0.05,0.95) length.lim <- c(28,72) ######################################## # Step 2 : call simulate function # # (Steps 3 to 6 occur within function) # ######################################## data.fin.df <- simulate.krill.call(No.ages.m,No.ages.f,M.f,M.mr,lambda,age0,No.ind,No.run,SD.m,SD.f,Avlen.m,Avlen.f, Mort.sd,Mort.m.sd,quant) # write results to csv file write.csv(x=data.fin.df, file = "Krill sex ratio simulation output 01.csv", na = "NA") # graph results print(xyplot(data.fin.df$Sex.r.q1+data.fin.df$Sex.r.q2+data.fin.df$Sex.r.q3+ rep(0.5,length(data.fin.df$length)) ~ data.fin.df$length, type=c("b","l","l","l"), lty=c(1,2,2,1),pch=c(1,NA,NA,NA), col=c(1,2,2,3), xlab="length (mm)", ylab="Sex Ratio (M/(M+F)", xlim=length.lim, ylim=c(-0.05,1.05))) print(xyplot(data.fin.df$Lfreq.m.q1+data.fin.df$Lfreq.m.q2+data.fin.df$Lfreq.m.q3+ data.fin.df$Lfreq.f.q1+data.fin.df$Lfreq.f.q2+data.fin.df$Lfreq.f.q3 ~ data.fin.df$length, type=rep(c("b","l","l"),2), lty=rep(c(1,2,2),2),pch=rep(c(1,NA,NA),2), col=c(1,1,1,4,4,4), xlab="length (mm)",ylab="Length Frequency (M & F)", xlim=length.lim, key = list(lines = list(col=c(1,1,1,4,4,4),lty=rep(c(1,2,2),2)), points=list(pch=rep(c(1,NA,NA),2),col=c(1,1,1,4,4,4)), text = list(lab = c("Males","lower 5%","upper 95%", "Females","lower 5%","upper 95%")),columns = 1, x=0.6,y=0.9))) print(xyplot(data.fin.df$Lfreq.t.q1+data.fin.df$Lfreq.t.q2+data.fin.df$Lfreq.t.q3 ~ data.fin.df$length, type=c("b","l","l"), lty=c(1,2,2),pch=c(1,NA,NA), col=c(1,2,2), xlab="length (mm)",ylab="Length Frequency (M+F)", xlim=length.lim)) ############################################################ # Functions to simulate krill sex ratio # ############################################################ # S.G. Candy # Vers 1.0 (12 Jan 2007) "simulate.krill.call" <- data.fin.df <- function(No.ages.m,No.ages.f,M.f,M.mr,lambda,age0,No.ind,No.run,SD.m,SD.f,Avlen.m,Avlen.f, Mort.sd,Mort.m.sd,quant) { # function to repeatedly call simulation using simulate.krill.sexr( ) and # summarise output # this version allows for accelerated male mortality # check correct number SD classes if (No.ages.m > length(SD.m)) print("not enough male SD classes") if (No.ages.f > length(SD.f)) print("not enough female SD classes") if (No.ages.m < length(SD.m)) print("too many male SD classes") if (No.ages.f < length(SD.f)) print("too many female SD classes") # graph survival function M.m <- M.f*M.mr P.surv.f <- c(1,exp(-M.f*(seq(2,No.ages.f,1)-1))) P.surv.m <- exp(-M.m*(seq(2,No.ages.m,1)-1))*((seq(2,No.ages.m,1)<=age0)+ (seq(2,No.ages.m,1)>age0)*exp(-(seq(2,No.ages.m,1)-age0)*lambda)) P.surv.m <- c(1,P.surv.m) if (No.ages.m < No.ages.f) P.surv.m <- c(P.surv.m,rep(0,(No.ages.f-No.ages.m))) if (No.ages.m > No.ages.f) P.surv.f <- c(P.surv.f,rep(0,(No.ages.m-No.ages.f))) No.ages <- max(c(No.ages.m,No.ages.f)) print(xyplot(P.surv.m+P.surv.f ~ seq(1,No.ages,1), type=c("l","l"), lty=c(1,2),col=c(1,2), xlab="age (yrs)", ylab="Survival Proportion", ylim=c(-0.05,1.05), key = list(lines = list(col=c(1,2),lty=c(1,2)), text = list(lab = c("Males","Females")),columns = 1, x=0.6,y=0.9))) # set up full set of length bin labels lev <- seq(20,70,1) nL <- length(lev)-1 lev.c <- character(length=nL+2) for (r in 1:nL) { lev.c[r+1] <- paste(paste("(",as.character(lev[r]),sep=""), paste(as.character(lev[r+1]),"]",sep=""),sep=",") } lev.c[1] <- "(0,20]" lev.c[nL+2] <- "(70,100]" length <- seq(19,70,1) Nlen <- length(length) Matn.len.f <- matrix(data=rep(0,No.run*Nlen), nrow=No.run, ncol=Nlen) Matn.len.m <- matrix(data=rep(0,No.run*Nlen), nrow=No.run, ncol=Nlen) Matn.len.sr <- matrix(data=rep(0,No.run*Nlen), nrow=No.run, ncol=Nlen) # call function No.run times and summarise results # if including stochastic variation in mortality if (Mort.sd < 0.0001) { for (i in 1:No.run) { data.df <- simulate.krill.sexr(No.ages.m,No.ages.f,M.f,M.m,lambda,age0,No.ind,SD.m,SD.f,Avlen.m,Avlen.f,lev,lev.c) Matn.len.m [i,] <- data.df$n.len.mf Matn.len.f [i,] <- data.df$n.len.ff Matn.len.sr [i,] <- data.df$sex.r }} if(Mort.sd > 0.0001) { # include variation in Mf for (i in 1:No.run) { M.fs <- M.f*exp(rnorm(n=1, mean=0, sd=Mort.sd)) M.ms <- M.mr*M.fs*exp(rnorm(n=1, mean=0, sd=Mort.m.sd)) data.df <- simulate.krill.sexr(No.ages.m,No.ages.f,M.fs,M.ms,lambda,age0,No.ind,SD.m,SD.f,Avlen.m,Avlen.f,lev,lev.c) Matn.len.m [i,] <- data.df$n.len.mf Matn.len.f [i,] <- data.df$n.len.ff Matn.len.sr [i,] <- data.df$sex.r }} # calculate quantiles (e.g. median, %5, 95%) and graph results Matq.len.f <- matrix(data=rep(0,3*Nlen), nrow=3, ncol=Nlen) Matq.len.m <- matrix(data=rep(0,3*Nlen), nrow=3, ncol=Nlen) Matq.len.t <- matrix(data=rep(0,3*Nlen), nrow=3, ncol=Nlen) Matq.len.sr <- matrix(data=rep(0,3*Nlen), nrow=3, ncol=Nlen) for (i in 1:Nlen) { Matq.len.f[,i] <- quantile(x=Matn.len.f [,i], probs = quant, na.rm = TRUE) Matq.len.m[,i] <- quantile(x=Matn.len.m [,i], probs = quant, na.rm = TRUE) Matq.len.sr[,i] <- quantile(x=Matn.len.sr [,i], probs = quant, na.rm = TRUE) Matq.len.t[,i] <- quantile(x=(Matn.len.m [,i]+Matn.len.f [,i]), probs = quant, na.rm = TRUE) } Sex.r.q1 <- Matq.len.sr[1,] Sex.r.q2 <- Matq.len.sr[2,] Sex.r.q3 <- Matq.len.sr[3,] Lfreq.m.q1 <- Matq.len.m[1,] Lfreq.m.q2 <- Matq.len.m[2,] Lfreq.m.q3 <- Matq.len.m[3,] Lfreq.f.q1 <- Matq.len.f[1,] Lfreq.f.q2 <- Matq.len.f[2,] Lfreq.f.q3 <- Matq.len.f[3,] Lfreq.t.q1 <- Matq.len.t[1,] Lfreq.t.q2 <- Matq.len.t[2,] Lfreq.t.q3 <- Matq.len.t[3,] data.fin.df <- data.frame(length,Sex.r.q1,Sex.r.q2,Sex.r.q3,Lfreq.m.q1,Lfreq.m.q2,Lfreq.m.q3, Lfreq.f.q1,Lfreq.f.q2,Lfreq.f.q3,Lfreq.t.q1,Lfreq.t.q2,Lfreq.t.q3) } "simulate.krill.sexr"<- krilldata.df <- function(No.ages.m,No.ages.f,M.f,M.m,lambda,age0,No.ind,SD.m,SD.f,Avlen.m,Avlen.f,lev,lev.c) { # function to run simulation on ages x No.ind nL <- length(lev)-1 # generate random normal deviates after thinning due to mortality # assume all year 1's survive so numbers are relative len.m <- rnorm(n=No.ind, mean=Avlen.m[1], sd=SD.m[1]) len.f <- rnorm(n=No.ind, mean=Avlen.f[1], sd=SD.f[1]) for (r in 2:No.ages.m){ No.surv.m <- No.ind*(exp(-M.m*(r-1))*(as.integer(r<=age0)+ as.integer(r>age0)*exp(-(r-age0)*lambda))) len.m <- c(len.m, rnorm(n=No.surv.m, mean=Avlen.m[r], sd=SD.m[r])) } for (r in 2:No.ages.f){ No.surv.f <- No.ind*exp(-M.f*(r-1)) len.f <- c(len.f,rnorm(n=No.surv.f, mean=Avlen.f[r], sd=SD.f[1])) } # aggregate into length bins len.m.f <- factor(x=cut(x=len.m, breaks=c(0,seq(20,70,1),100))) len.f.f <- factor(x=cut(x=len.f, breaks=c(0,seq(20,70,1),100))) n.len.m <- as.vector(tapply(X=rep(1,length(len.m)), INDEX=len.m.f, FUN=sum)) n.len.f <- as.vector(tapply(X=rep(1,length(len.f)), INDEX=len.f.f, FUN=sum)) m.ind <- (seq(1,(nL+2)))[lev.c %in% (levels(len.m.f))] n.len.mf <- rep(0,(nL+2)) n.len.mf[m.ind] <- n.len.m f.ind <- (seq(1,(nL+2)))[lev.c %in% (levels(len.f.f))] n.len.ff <- rep(0,(nL+2)) n.len.ff[f.ind] <- n.len.f n.tot <- n.len.mf + n.len.ff n.tot[n.tot==0] <- NA sex.r <- rep(1,length(n.tot)) sex.r[!is.na(n.tot)] <- n.len.mf/n.tot sex.r[is.na(n.tot)] <- NA length <- seq(19,70,1) krilldata.df <- data.frame(sex.r,length,n.len.mf,n.len.ff) }