#!/usr/bin/env python # coding: utf-8 # In[1]: get_ipython().run_line_magic('matplotlib', 'inline') import pymc3 as pm import numpy as np import pandas as pd import theano.tensor as tt from theano import shared, config from theano.tensor.shared_randomstreams import RandomStreams srng = RandomStreams(seed=20090425) # In[23]: config.exception_verbosity='high' # In[3]: data15 = pd.read_csv('mortalities.15km.dword.1996.2013.csv') data30 = pd.read_csv('mortalities.30km.dword.1996.2013.csv') # In[4]: # FL regions regions = ['ATL', 'USJ', 'NW', 'SW'] ATL, USJ, NW, SW = range(4) # Winter habitat quality habitats = ['Low', 'Medium', 'High'] LOW, MED, HIGH = range(3) # Causes of death causes = ['watercraft', 'wcs', 'debris', 'cold', 'redtide', 'other'] long_causes = causes + ['undetermined'] WATERCRAFT, WCS, DEBRIS, COLD, REDTIED, OTHER, UNDETERMINED = range(7) # Age classes classes = ['Calves', 'Subadults', 'Adults'] STARTYEAR, ENDYEAR = 1996, 2013 # Unusual mortality event red tide years tide_years = (1996, 2002, 2003, 2005, 2006, 2012, 2013) # Moderately unusual red tide years tide_years0 = (2002, 2003, 2005, 2006, 2012) # Intense red tide years tide_years_intense = (1996, 2013) # Atlantic unusual mortality event red tide years atl_tide_years = (2007, 2008) # Severe cold years severe_years = {ATL: range(2010, 2011), USJ: range(2010, 2011), NW: (2003, 2010), SW: range(2010, 2011)} # Cold years cold_years = {ATL: (1996, 2001, 2003), USJ: (1996, 2001, 2009, 2011), NW: (1996, 1998, 2001, 2011), SW: (1996, 2001)} # Uninformative priors for unknown proportions # In[5]: ncauses = len(causes), len(causes)-1, len(causes), len(causes) prior0 = [np.ones(k) for k in ncauses] # From Langtimm et al. 2004 (mortality ratios from USJ) # In[6]: calf1_mort_ratio = 0.190/0.031 calf2_mort_ratio = 0.085/0.031 # Adult baseline mortalities and se from Langtimm et al. unpublished (2015). # # Changed this to a tuple, so that you can just index with the appropriate constant # In[7]: adult_mort_base = 0.027417, 0.020739, 0.022355, 0.023073 adult_mort_base[SW] # Generate survival rates and variances on logit scale # In[8]: calf1_mort_base = np.array(adult_mort_base) * calf1_mort_ratio calf2_mort_base = np.array(adult_mort_base) * calf2_mort_ratio calf_mort_base = 1 - np.sqrt((1 - calf1_mort_base)*(1 - calf2_mort_base)) wcs_protect_year = 2001 # In[9]: calf_mort_base # In[10]: analysis_subset = (data30[(data30.year >= STARTYEAR) & (data30.year <= ENDYEAR) & (data30.age == 'Calves')] .replace({'area': dict(zip(regions, range(4))), 'habitat': dict(zip(habitats, range(3)))})) # In[11]: analysis_subset.head() # In[24]: def mortality_model(data, prior=prior0, verbose=0): # Jeff: simplifying by getting rid of ATL red tide, I'm not really planning to use that anymore. with pm.Model(verbose=verbose) as model: # Red tide effect factors a_mort = pm.Uniform('a_mort', lower=0, upper=1-calf_mort_base[SW]) a = a_mort / calf_mort_base[SW] b_mort = pm.Uniform('b_mort', lower=0, upper=1-calf_mort_base[SW]) b = b_mort / calf_mort_base[SW] # Winter effect factors for each habitat type sev_mort = pm.Uniform('sev_mort', lower=0, upper=1-max(calf_mort_base), shape=(3,4)) cold_mort_upper = np.resize(1-max(calf_mort_base), (3,4)) cold_mort_upper[HIGH] = 0.0001 cold_mort = pm.Uniform('cold_mort', lower=0, upper=cold_mort_upper, shape=(3,4)) norm_mort_upper = np.resize(1-max(calf_mort_base), (3,4)) norm_mort_upper[HIGH] = 0.0001 norm_mort = pm.Uniform('norm_mort', lower=0, upper=norm_mort_upper, shape=(3,4)) sev = sev_mort / np.resize(calf_mort_base, (3,4)) cold = cold_mort / np.resize(calf_mort_base, (3,4)) norm = norm_mort / np.resize(calf_mort_base, (3,4)) p = [pm.Dirichlet('p_{}'.format(a), t) for t,a in zip(prior, regions)] π = [pm.Dirichlet('π_{}'.format(a), np.ones(k)) for k,a in zip(ncauses, regions)] for (area, year, habitat), data in data.groupby(['area', 'year', 'habitat']): undetermined = data.undetermined.values Z = (data[causes].drop('redtide', axis=1).values if area==USJ else data[causes].values).squeeze() # Only moderately unusual red tide (no cold effects) if ((year in tide_years0) and (area==SW) and (habitat==HIGH)): # theta varies from pi based on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + a) / (1+a) η = p[area] # Moderately unusual red tide and normal cold elif ((year in tide_years0) and (area==SW)): no = norm[habitat, area] # theta varies from pi based on winter severity and on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + a) θ = tt.set_subtensor(θ[3], θ[3] + no) / (1 + a + no) # eta varies from p based on winter severity η = tt.set_subtensor(p[area][3], p[area][3] + no) / (1 + no) # Intense unusual red tide and cold year elif ((year in tide_years_intense) and (area==SW) and (year in cold_years[area])): co = cold[habitat, area] # theta varies from pi based on winter severity and on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + b) θ = tt.set_subtensor(θ[3], θ[3] + co) / (1 + b + co) # eta varies from p based on winter severity η = tt.set_subtensor(p[area][3], p[area][3] + co) / (1 + co) # Intense unusual red tide and severe cold year elif ((year in tide_years_intense) and (area==SW) and (year in severe_years[area])): sv = sev[habitat, area] # theta varies from pi based on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + b) θ = tt.set_subtensor(θ[3], θ[3] + sv) / (1 + b + sv) # eta varies from p based on winter severity η = tt.set_subtensor(p[area][3], p[area][3] + sv) / (1 + sv) # Intense unusual red tide only elif ((year in tide_years_intense) and (area==SW) and (habitat==HIGH)): # theta varies from pi based on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + b) / (1+b) η = p[area] # Intense unusual red tide and normal cold elif (year in tide_years_intense) and (area==SW): no = norm[habitat, area] # theta varies from pi based on winter severity and on red tide UME θ = tt.set_subtensor(π[area][4], π[area][4] + b) θ = tt.set_subtensor(θ[3], θ[3] + no) / (1 + b + no) # eta varies from p based on winter severity η = tt.set_subtensor(p[area][3], p[area][3] + no) / (1 + no) # Severe cold year elif (year in severe_years[area]): sv = sev[habitat, area] # theta varies from pi based on winter severity θ = tt.set_subtensor(π[area][3], π[area][3] + sv) / (1 + sv) # eta varies from p based on winter severity η = tt.set_subtensor(p[area][3], p[area][3] + sv) / (1 + sv) # Cold year, low or medium habitat quality elif (year in cold_years[area]) and (habitat