#This routine calculates some basic risk metrics using the DennisHolmes method (aka Holmes and Fagan 2002) #Confidence intervals on mu, Probability of extinction with 95% confidence intervals, #Probability of decline with 95% confidence intervals, Probability of observing extinction, #Probability of observing a given decline # #mu: estimate of parameter mu #sigma2: estimate of parameter sigma2_p=sigma2_slp #sigma2np: estimate of paramter sigma2_np #nq: current population size; absolute for extinction calcs or relative for decline calcs #ne: forecast population size; absolute for extinction calcs or relative for decline calcs #n: number of counts #L: filter length #years: number of forecast years; used in probability of decline and extinction calcs #probofext: F means do not do prob of extinction calculations #probofdecline: F means do not do prob of decline calculations # #Example 1: To calculate the probability of 90% decline in 25 years for a population # riskcalcs(mu=-.05, sigma2=.05, nq=10, ne=1, n=20, L=4, years=25, probofdecline=T) #Example 2: To calculate the probability of extinction in 25 years for a population of size 200 # riskcalcs(mu=-.05, sigma2=.05, nq=200, ne=1, n=20, L=4, years=25, probofext=T) #Example 3: To calculate the probability of OBSERVING extinction in 25 years for a population of size 200 # riskcalcs(mu=-.05, sigma2=.05, sigma2np=.2, nq=200, ne=1, n=20, L=4, years=25, probofext=T) riskcalcs_function(mu=NA, sigma2=NA, sigma2np=0, nq=NA, ne=NA, n=NA, L=NA, years=1, probofext=F, probofdecline=F) { if(is.na(mu)){print("Error denniscalcs; need mu"); break} if(is.na(sigma2)){print("Error denniscalcs; need sigma2"); break} if(is.na(nq)){print("Error denniscalcs; need nq"); break} if(is.na(ne)){print("Error denniscalcs; need ne"); break} if(is.na(n) | is.na(L)){print("Error denniscalcs; need n and L"); break} if(n<10) {print("Error denniscalcs; n must be greater than 10"); break} if(sigma2==0) sigma2_.000001 #will choke on 0 tq_n-L+1 if(n==10) dfRslp_1.42 else dfslp_.333+.212*n-.387*L #from Holmes and Fagan 2002 q _ tq-1; d _ .000001 cumprobext_NA cumprobextwithobserror_NA cumprobdec_NA cumprobdecwithobserror_NA cumprobext95confints_NA cumprobdec95confints_NA mu95confints_NA #Get confidence intervals on lambda muup95 _ ( mu - qt(.025,dfslp)*sqrt(sigma2*(1/q)) ) mulow95 _ ( mu + qt(.025,dfslp)*sqrt(sigma2*(1/q)) ) mu95confints _ c(muup95,mulow95) #Get extinction probabilities on true population size if(probofext) { tmp_2*log(nq/ne)*abs(mu)/sigma2 tmp_tmp/700 #this is a trick to get rid of NAs when tmp is big if(tmp<1)tmp_1 G _ pnorm( (-log(nq/ne) + abs(mu)*years)/sqrt(sigma2*years)) + (exp((1/tmp)*2*log(nq/ne)*abs(mu)/sigma2)* pnorm( (-log(nq/ne) - abs(mu)*years)/sqrt(sigma2*years))^(1/tmp))^tmp; if(mu < 0) probext _1 else probext _ (ne/nq)^(2*mu/(sigma2*(q-1)/q)) cumprobext _ G*probext; } #Get probability true population is nq/ne lower at time years if (probofdecline) #Prob that stock is nq/ne lower at year t { if(sigma2np < 0) sigma2np_0 G _ 1-pnorm( (log(nq/ne) + mu*years)/sqrt(sigma2*years)); cumprobdec _ G; } #Get confidence intervals on prob of decline or extinction via monte carlo simulation following dennis 1991 if (probofext | probofdecline) { tmpreps_1000 murnds_rt(tmpreps,dfslp)*sqrt(sigma2/(n-L))+mu muhigh_qt(.975,dfslp)*sqrt(sigma2/(n-L))+mu mulow_qt(.025,dfslp)*sqrt(sigma2/(n-L))+mu s2slprnds_sigma2*dfslp/rchisq(tmpreps,dfslp) s2low_sigma2*dfslp/qchisq(.975,dfslp) s2high_sigma2*dfslp/qchisq(.025,dfslp) okparams_(murnds<=muhigh) & (murnds>=mulow) & (s2slprnds<=s2high) & (s2slprnds>=s2low) murnds_murnds[okparams] s2slprnds_s2slprnds[okparams] if(probofext) { cumprobext95confints_{} for (i in 1:length(years)) { tmp_2*log(nq/ne)*abs(murnds)/s2slprnds tmp_tmp/700 #this is a trick to get rid of NAs when tmp is big tmp[tmp<1]_1 G _ pnorm( (-log(nq/ne) + abs(murnds)*years[i])/sqrt(s2slprnds*years[i])) + (exp((1/tmp)*2*log(nq/ne)*abs(murnds)/s2slprnds)* pnorm( (-log(nq/ne) - abs(murnds)*years[i])/sqrt(s2slprnds*years[i]))^(1/tmp))^tmp; probext _ (ne/nq)^(2*murnds/(s2slprnds*(q-1)/q)) probext[murnds<0]_1 cumprobextrnds _ G*probext; upper95 _ quantile(cumprobextrnds,.975); low95 _ quantile(cumprobextrnds,.025); cumprobext95confints_c(cumprobext95confints,quantile(cumprobextrnds,c(.975,.025))) } }#if probofext=T if(probofdecline) { cumprobdec95confints_{} for (i in 1:length(years)) { cumprobdecrnds _ 1-pnorm( (log(nq/ne) + murnds*years[i])/sqrt(s2slprnds*years[i])); cumprobdec95confints_c(cumprobdec95confints,quantile(cumprobdecrnds,c(.975,.025))) } }#if probofdecline=T } #probconfints #Probability of observing the stock is below threshold with time tau given that observations have error if (probofext==T & sigma2np!=0) { prob_0 reps_500 tau_1/1 #basically the step length projtime_years/tau for(nn in 1:reps) { N_rep(0,projtime+1) N[1]_100 ep_exp(rnorm(projtime+1,mu*tau,sqrt(sigma2*tau))) if(sigma2np<=0) enp_1 else enp_exp(rnorm(projtime+1,0,sqrt(sigma2np))) for(year in 2:(projtime+1)) N[year]_N[year-1]*ep[year] N_enp*N if(min(N[seq(1/tau,projtime,1/tau)])<=(ne/nq)*N[1]) prob_prob+1 } cumprobextwithobserror_prob/reps } #Prob of observing stock below threshold at time tau given that observations are corrupted if (probofdecline==T & sigma2np!=0) #Prob that stock is nq/ne lower at year t { reps_500 projtime_years prob_0 for(nn in 1:reps) { N_rep(0,projtime+1) N[1]_100 ep_exp(rnorm(projtime+1,mu,sqrt(sigma2))) if(sigma2np<=0) enp_1 else enp_exp(rnorm(projtime+1,0,sqrt(sigma2np))) for(year in 2:(projtime+1)) N[year]_N[year-1]*ep[year] N_enp*N if(N[projtime+1]<=(ne/nq)*N[1]) prob_prob+1 } cumprobdecwithobserror_prob/reps } output _ list(mu95confints=mu95confints, cumprobext=cumprobext, cumprobext95confints=cumprobext95confints, cumprobdec=cumprobdec, cumprobdec95confints=cumprobdec95confints, cumprobextwithobserror=cumprobextwithobserror, cumprobdecwithobserror=cumprobdecwithobserror ) return(output) }