data = read.csv(file.choose()) #use this to read in "TKW_presence.csv" ids = data$id #pull out IDs and age-sex classes asc = data$asc #vector to hold results of simulations seals_eaten = rep(NA,10000) #prey preferences from literature (see text) seal_p = 0.69 steller_p = 0.12 porp_p = 0.13 dall_p = 0.06 for(i in 1:10000){ #generate average caloric content of prey items, using ranges from literature (see text) steller_content = 2500*runif(1,420,725) dall_content = 4730*runif(1,170,200) porp_content = 4730*runif(1,55,70) seal_content = 3500*runif(1,50,100) seal_portion = (seal_content*seal_p)/(seal_content*seal_p + steller_content*steller_p + porp_content*porp_p + dall_content*dall_p) #calculate the portion of energetic need met by seals kw_mass = rep(NA,length(ids)) #vector to hold masses kw_mass[asc == "AM"] = runif(sum(asc=="AM"),4000,8000) #randomly generate masses within specified range (see text) kw_mass[asc == "AF"] = runif(sum(asc=="AF"),3000,4000) kw_mass[asc == "IM"] = runif(sum(asc=="IM"),1500,2000) kw_energy = ((19.65*kw_mass^0.756)*0.86)*24 #calculate energy needs in kcal/d kw_seal_kcal = kw_energy*seal_portion #how many kcal/d of seal are needed seal_available = seal_content*0.85 #factor in assimilation value seals_daily = kw_seal_kcal/seal_available #change to number of seals seals_eaten[i] = sum(seals_daily*data$days_in_salish_sea) #total seals eaten by all whales } mean_estimate = mean(seals_eaten) confidence_interval = quantile(seals_eaten, c(0.025,0.975))