Supplement 1 - All Python/STAN Code, Including Bayesian Hierarchical Model, Used in the Analysis. Nathan P. Lemoine (1,5), Jessica Shue (2), Brittany Verrico (3), David Erickson (4), W. John Kress (4), John D. Parker (2) 1 Department of Biological Sciences, Florida International University North Miami, FL 33181 2 Smithsonian Environmental Research Center, Smithsonian Institution Edgewater, MD 33181 3 Washington and Jefferson College, Washington, PA 15301 4 National Museum of Natural History, Smithsonian Institution Washington DC, 20560 5 Current Affiliation: Department of Biology, Colorado State University Fort Collins, CO 80523 Correspondence: Nathan P. Lemoine Department of Biology Colorado State University Room E106 Anatomy/Zoology Building Fort Collins, CO 80523 phone: 305-919-4109, email: lemoine.nathan@gmail.com #### CODE STARTS HERE import pandas as pd import pylab as py import numpy as np import ete2 as ete import scipy.spatial.distance as dst import scipy.stats as st import pystan # load in the data traits = pd.read_csv('/Users/Nate/Documents/FIU/Research/Invasion_TraitPhylo/Data/plantTraits.csv') surveyData = pd.read_csv('/Users/Nate/Documents/FIU/Research/Invasion_TraitPhylo/Data/master_sp_data.csv') envData = pd.read_csv('/Users/Nate/Documents/FIU/Research/Invasion_TraitPhylo/Data/master_env_data.csv') SERCphylo = ete.Tree('/Users/Nate/Documents/FIU/Research/SERC_Phylo/SERC_Nov1-2013.newick.tre') #### TRAIT CLEANUP #### # put an underscore in trait species traits['species'] = traits['species'].map(lambda x: x.replace(' ', '_')) # pull out the relevant traits and only keep complete cases traits = traits[['species', 'introduced', 'woody', 'SLA', 'seedMass', 'toughness']] traits = traits.dropna() # standardize the traits traitsSTD = traits traitsSTD[['SLA', 'seedMass', 'toughness']] = traitsSTD[['SLA', 'seedMass', 'toughness']].apply(lambda x: (x - x.mean())/x.std()) #### ENVIRONMENTAL CLEANUP #### # summarize by plot envPlot = envData[['plot', 'avgvwc', 'trans', 'litter', 'age']].groupby('plot', as_index = False).mean() envPlot['light_scale'] = (envPlot['trans'] - envPlot['trans'].mean()) / envPlot['trans'].std() envPlot['water_scale'] = (envPlot['avgvwc'] - envPlot['avgvwc'].mean()) / envPlot['avgvwc'].std() envPlot['litter_scale'] = (envPlot['litter'] - envPlot['litter'].mean()) / envPlot['litter'].std() envPlot['age_scale'] = (envPlot['age'] - envPlot['age'].mean()) / envPlot['age'].std() #### COVER DATA CLEANUP #### # insert underscore in cover species surveyData['species'] = surveyData['species'].map(lambda x: x.replace(' ', '_')) coverData = surveyData[surveyData['species'].isin(traitsSTD['species'])] # fill in the NA's in diagonal with filler (otherwise they get excluded later on coverData['diagonal'].fillna(value = 'N', inplace = True) # put into wide format quadSpMat = coverData.pivot_table(values = 'cover', rows = ['plot', 'diagonal','quadrat'], cols = 'species', aggfunc = 'mean') # fill the NA's with 0 quadSpMat.fillna(value = 0, inplace = True) # sum over plots standSpMat = quadSpMat.groupby(level = 0).sum() # melt back into long format standSum = standSpMat.stack() # calculate relative cover standRelCover = standSum.groupby(level = 0).apply(lambda x: x/x.sum()) standRelCover = standRelCover.reset_index() standRelCover.rename(columns = {0 : 'relCover'}, inplace = True) # define the actual function that calculates MPD for the species def MPDcalc(sp, tempSpecies, phylo): # get the name of the species currentSp = sp # drop that species from the data frame templist = tempSpecies[tempSpecies['species'] != currentSp] # calculate MPD MPD = templist['species'].apply(lambda a: phylo.get_distance(currentSp, a)) # weight MPD by relative abundance MPD = MPD * tempSpecies['relCover'] return MPD.sum() # define a function that calculates MPD and MTD for every species def distCalc(group): # remove species not in the phylogeny tempSpecies = group[group['species'].isin(SERCphylo.get_leaf_names())] # prune the phylogeny to be site specific by first making a copy of the SERC phylo and then pruning the copy tempPhylo = SERCphylo.copy() tempPhylo.prune(tempSpecies['species'], preserve_branch_length = True) # for every species, calculate the mean pairwise distance tempSpecies['MPD'] = tempSpecies['species'].apply(lambda x: MPDcalc(x, tempSpecies, tempPhylo)) # get traits for species present at the site and put them in the correct order so that relative abundances line up correctly! tempTraits = traitsSTD[traitsSTD['species'].isin(tempSpecies['species'])] tempTraits.sort('species', inplace = True) # calculate pairwise distances distances = dst.squareform( dst.pdist(tempTraits[['SLA', 'seedMass', 'toughness']], 'euclidean') ) # calculate trait-weighted MPD (this is akin to multiplying each distance by relative abundance and summing since the weights add to 1) tempSpecies['MTD'] = distances.dot(tempSpecies['relCover']) return tempSpecies # run the function fullData = standRelCover.groupby('plot', as_index = False).apply(distCalc) # next, merge the species with the introduced/native status fullData = pd.merge(fullData, traitsSTD[['species', 'introduced']]) # compare the MPD of native and introduced species intGroup = fullData.groupby('introduced') intGroup['MPD'].agg({'mean' : np.mean, 'sd' : lambda x: x.std()}) # compare the MTD intGroup['MTD'].agg({'mean' : np.mean, 'sd' : lambda x: x.std()}) # turn the plot ID into a numeric factor (the environmental plot data is already in order) fullData['plotID'] = pd.factorize(fullData['plot'])[0] + 1 # turn Presence/Absence into a 1,0 column fullData['PA'] = (fullData['relCover'] > 0).astype('int') # plot histogram of relative cover values py.hist(fullData['relCover'], bins = 50, normed = True) py.show() # plot histogram of relative cover values (only those > 0 ) py.hist(fullData[fullData['relCover'] > 0]['relCover'], bins = 50, normed = True, alpha = 0.3) # fit a gamma distribution param = st.gamma.fit(fullData[fullData['relCover'] > 0]['relCover']) x = np.linspace(0.0004, 1, 100) pdf_fitted = st.gamma.pdf(x, param[0], param[1], param[2]) py.plot(x, pdf_fitted, 'r--', label = 'Gamma Fit') py.xlabel('Relative Cover') py.ylabel('Density') py.legend() py.savefig('gamma_hist.pdf', bbox_inches = 'tight') py.show() # scale MPD and MTD fullData['MPD_scale'] = (fullData['MPD'] - fullData['MPD'].mean()) / fullData['MPD'].std() fullData['MTD_scale'] = (fullData['MTD'] - fullData['MTD'].mean()) / fullData['MTD'].std() fullData.sort(['plot', 'species'], inplace = True) # set up the STAN model hurdleMod = """ data{ int N; // number of observations int J; // number of plots vector[N] y; // relative cover vector[N] MPD; // MPD values vector[N] MTD; // MTD values real Origin[N]; // Origin int plot[N]; // Plot ID vector[J] light; // Plot light vector[J] vwc; // Plot VWC vector[J] litter; // Plot litter vector[J] age; // Plot age matrix[6,6] R; // Prior on Taus } parameters{ vector[6] beta[J]; // presence/absence coefficients vector[6] kappa[J]; // relative cover coefficients real shape; // shape parameter // PRESENCE/ABSENCE GROUP-LEVEL COEFFICIENTS real int_mu; real mpd_mu; real mtd_mu; real origin_mu; real mpdOrigin_mu; real mtdOrigin_mu; real int_light; real mpd_light; real mtd_light; real origin_light; real mpdOrigin_light; real mtdOrigin_light; real int_vwc; real mpd_vwc; real mtd_vwc; real origin_vwc; real mpdOrigin_vwc; real mtdOrigin_vwc; real int_litter; real mpd_litter; real mtd_litter; real origin_litter; real mpdOrigin_litter; real mtdOrigin_litter; real int_age; real mpd_age; real mtd_age; real origin_age; real mpdOrigin_age; real mtdOrigin_age; // RELATIVE COVER GROUP LEVEL COEFFICIENTS real int_mu_RC; real mpd_mu_RC; real mtd_mu_RC; real origin_mu_RC; real mpdOrigin_mu_RC; real mtdOrigin_mu_RC; real int_light_RC; real mpd_light_RC; real mtd_light_RC; real origin_light_RC; real mpdOrigin_light_RC; real mtdOrigin_light_RC; real int_vwc_RC; real mpd_vwc_RC; real mtd_vwc_RC; real origin_vwc_RC; real mpdOrigin_vwc_RC; real mtdOrigin_vwc_RC; real int_litter_RC; real mpd_litter_RC; real mtd_litter_RC; real origin_litter_RC; real mpdOrigin_litter_RC; real mtdOrigin_litter_RC; real int_age_RC; real mpd_age_RC; real mtd_age_RC; real origin_age_RC; real mpdOrigin_age_RC; real mtdOrigin_age_RC; cov_matrix[6] Tau_Beta; cov_matrix[6] Tau_Kappa; } transformed parameters{ vector[N] theta; // logit predictor vector[N] pi; // conversion of logit to probability vector[N] eta; // gamma predictor vector[N] yhat; // exponentiate gamma link vector[N] rate; // rate parameter for the gamma distritbution vector[6] betaHat[J]; vector[6] kappaHat[J]; // DEFINE THE HURDLE MODEL for(n in 1:N){ // PRESENCE ABSENCE theta[n] <- beta[plot[n],1] + beta[plot[n],2]*MPD[n] + beta[plot[n],3]*MTD[n] + beta[plot[n],4]*Origin[n] + beta[plot[n],5]*MPD[n]*Origin[n] + beta[plot[n],6]*MTD[n]*Origin[n]; pi[n] <- inv_logit(theta[n]); // RELATIVE COVER eta[n] <- kappa[plot[n],1] + kappa[plot[n],2]*MPD[n] + kappa[plot[n],3]*MTD[n] + kappa[plot[n],4]*Origin[n] + kappa[plot[n],5]*MPD[n]*Origin[n] + kappa[plot[n],6]*MTD[n]*Origin[n]; yhat[n] <- exp(eta[n]); rate[n] <- shape / yhat[n]; } // GROUP-LEVEL MODEL WHERE RANDOM EFFECTS ARE FUNCTIONS OF GROUP-LEVEL PREDICTORS for(j in 1:J){ betaHat[j, 1] <- int_mu + int_light*light[j] + int_vwc*vwc[j] + int_litter*litter[j] + int_age*age[j]; betaHat[j, 2] <- mpd_mu + mpd_light*light[j] + mpd_vwc*vwc[j] + mpd_litter*litter[j] + mpd_age*age[j]; betaHat[j, 3] <- mtd_mu + mtd_light*light[j] + mtd_vwc*vwc[j] + mtd_litter*litter[j] + mtd_age*age[j]; betaHat[j, 4] <- origin_mu + origin_light*light[j] + origin_vwc*vwc[j] + origin_litter*litter[j] + origin_age*age[j]; betaHat[j, 5] <- mpdOrigin_mu + mpdOrigin_light*light[j] + mpdOrigin_vwc*vwc[j] + mpdOrigin_litter*litter[j] + mpdOrigin_age*age[j]; betaHat[j, 6] <- mtdOrigin_mu + mtdOrigin_light*light[j] + mtdOrigin_vwc*vwc[j] + mtdOrigin_litter*litter[j] + mtdOrigin_age*age[j]; kappaHat[j, 1] <- int_mu_RC + int_light_RC*light[j] + int_vwc_RC*vwc[j] + int_litter_RC*litter[j] + int_age_RC*age[j]; kappaHat[j, 2] <- mpd_mu_RC + mpd_light_RC*light[j] + mpd_vwc_RC*vwc[j] + mpd_litter_RC*litter[j] + mpd_age_RC*age[j]; kappaHat[j, 3] <- mtd_mu_RC + mtd_light_RC*light[j] + mtd_vwc_RC*vwc[j] + mtd_litter_RC*litter[j] + mtd_age_RC*age[j]; kappaHat[j, 4] <- origin_mu_RC + origin_light_RC*light[j] + origin_vwc_RC*vwc[j] + origin_litter_RC*litter[j] + origin_age_RC*age[j]; kappaHat[j, 5] <- mpdOrigin_mu_RC + mpdOrigin_light_RC*light[j] + mpdOrigin_vwc_RC*vwc[j] + mpdOrigin_litter_RC*litter[j] + mpdOrigin_age_RC*age[j]; kappaHat[j, 6] <- mtdOrigin_mu_RC + mtdOrigin_light_RC*light[j] + mtdOrigin_vwc_RC*vwc[j] + mtdOrigin_litter_RC*litter[j] + mtdOrigin_age_RC*age[j]; } } model{ // ACTUAL HURDLE MODEL for(n in 1:N){ increment_log_prob( if_else(y[n] == 0, log1m(pi[n]), log(pi[n]) + gamma_log(y[n], shape, rate[n])) ); } // MULTI-NORMAL DISTRIBUTION FOR RANDOM EFFECTS for(j in 1:J){ beta[j] ~ multi_normal_prec(betaHat[j], Tau_Beta); kappa[j] ~ multi_normal_prec(kappaHat[j], Tau_Kappa); } // PRIORS - PUT NONINFORMATIVE PRIORS ON THE GROUP-LEVEL PARAMETERS, WISHART DISTRIBUTIONS FOR THE PRECISIONS // USE NORMAL(0, 2) TO SPEED CONVERGENCE AND PREVENT THE CHAINS FROM WANDERING TOO FAR AFIELD (ALL PREDICTORS STANDARDIZED) int_mu ~ normal(0, 2); mpd_mu ~ normal(0, 2); mtd_mu ~ normal(0, 2); origin_mu ~ normal(0, 2); mpdOrigin_mu ~ normal(0, 2); mtdOrigin_mu ~ normal(0, 2); int_light ~ normal(0, 2); mpd_light ~ normal(0, 2); mtd_light ~ normal(0, 2); origin_light ~ normal(0, 2); mpdOrigin_light ~ normal(0, 2); mtdOrigin_light ~ normal(0, 2); int_vwc ~ normal(0, 2); mpd_vwc ~ normal(0, 2); mtd_vwc ~ normal(0, 2); origin_vwc ~ normal(0, 2); mpdOrigin_vwc ~ normal(0, 2); mtdOrigin_vwc ~ normal(0, 2); int_litter ~ normal(0, 2); mpd_litter ~ normal(0, 2); mtd_litter ~ normal(0, 2); origin_litter ~ normal(0, 2); mpdOrigin_litter ~ normal(0, 2); mtdOrigin_litter ~ normal(0, 2); int_age ~ normal(0, 2); mpd_age ~ normal(0, 2); mtd_age ~ normal(0, 2); origin_age ~ normal(0, 2); mpdOrigin_age ~ normal(0, 2); mtdOrigin_age ~ normal(0, 2); int_mu_RC ~ normal(0, 2); mpd_mu_RC ~ normal(0, 2); mtd_mu_RC ~ normal(0, 2); origin_mu_RC ~ normal(0, 2); mpdOrigin_mu_RC ~ normal(0, 2); mtdOrigin_mu_RC ~ normal(0, 2); int_light_RC ~ normal(0, 2); mpd_light_RC ~ normal(0, 2); mtd_light_RC ~ normal(0, 2); origin_light_RC ~ normal(0, 2); mpdOrigin_light_RC ~ normal(0, 2); mtdOrigin_light_RC ~ normal(0, 2); int_vwc_RC ~ normal(0, 2); mpd_vwc_RC ~ normal(0, 2); mtd_vwc_RC ~ normal(0, 2); origin_vwc_RC ~ normal(0, 2); mpdOrigin_vwc_RC ~ normal(0, 2); mtdOrigin_vwc_RC ~ normal(0, 2); int_litter_RC ~ normal(0, 2); mpd_litter_RC ~ normal(0, 2); mtd_litter_RC ~ normal(0, 2); origin_litter_RC ~ normal(0, 2); mpdOrigin_litter_RC ~ normal(0, 2); mtdOrigin_litter_RC ~ normal(0, 2); int_age_RC ~ normal(0, 2); mpd_age_RC ~ normal(0, 2); mtd_age_RC ~ normal(0, 2); origin_age_RC ~ normal(0, 2); mpdOrigin_age_RC ~ normal(0, 2); mtdOrigin_age_RC ~ normal(0, 2); Tau_Beta ~ wishart(7, R); Tau_Kappa ~ wishart(7, R); } """ data = {'N' : len(fullData), 'J' : fullData['plotID'].max(), 'plot' : fullData['plotID'], 'y' : fullData['relCover'], 'MPD' : fullData['MPD_scale'], 'MTD' : fullData['MTD_scale'], 'Origin' : fullData['introduced'], 'light' : envPlot['light_scale'], 'vwc' : envPlot['water_scale'], 'litter' : envPlot['litter_scale'], 'age' : envPlot['age_scale'], 'R' : np.identity(6)} fit = pystan.stan(model_code = hurdleMod, data = data, chains = 4, iter = 1000) ## check convergance of ALL parameters # first, intercept fit.plot(['int_mu', 'int_light', 'int_vwc', 'int_litter', 'int_age']) fit.plot(['int_mu_RC', 'int_light_RC', 'int_vwc_RC', 'int_litter_RC', 'int_age_RC']) # next, MPD fit.plot(['mpd_mu', 'mpd_light', 'mpd_vwc', 'mpd_litter', 'mpd_age']) fit.plot(['mpd_mu_RC', 'mpd_light_RC', 'mpd_vwc_RC', 'mpd_litter_RC', 'mpd_age_RC']) # next, MTD fit.plot(['mtd_mu', 'mtd_light', 'mtd_vwc', 'mtd_litter', 'mtd_age']) fit.plot(['mtd_mu_RC', 'mtd_light_RC', 'mtd_vwc_RC', 'mtd_litter_RC', 'mtd_age_RC']) # next, Origin fit.plot(['origin_mu', 'origin_light', 'origin_vwc', 'origin_litter', 'origin_age']) fit.plot(['origin_mu_RC', 'origin_light_RC', 'origin_vwc_RC', 'origin_litter_RC', 'origin_age_RC']) # next, MPD:Origin fit.plot(['mpdOrigin_mu', 'mpdOrigin_light', 'mpdOrigin_vwc', 'mpdOrigin_litter', 'mpdOrigin_age']) fit.plot(['mpdOrigin_mu_RC', 'mpdOrigin_light_RC', 'mpdOrigin_vwc_RC', 'mpdOrigin_litter_RC', 'mpdOrigin_age_RC']) # next, MTD:Origin fit.plot(['mtdOrigin_mu', 'mtdOrigin_light', 'mtdOrigin_vwc', 'mtdOrigin_litter', 'mtdOrigin_age']) fit.plot(['mtdOrigin_mu_RC', 'mtdOrigin_light_RC', 'mtdOrigin_vwc_RC', 'mtdOrigin_litter_RC', 'mtdOrigin_age_RC']) py.show() # next, the betas and kappas # first make a function that does the density plotting for each dataframe def DensPlot(data): fig, axes = py.subplots(nrows = 5, ncols = 5) # make a subplots object for i in range(25): ix = np.unravel_index(i, axes.shape) # get the subplot location (i.e. (0,0)) from the integer index data[data.columns[i]].plot(kind = 'kde', ax = axes[ix]) py.show() # next make a function that does the trace plotting for each dataframe def TracePlot(data): fig, axes = py.subplots(nrows = 5, ncols = 5) # make a subplots object for i in range(25): ix = np.unravel_index(i, axes.shape) # get the subplot location (i.e. (0,0)) from the integer index data[data.columns[i]].plot(ax = axes[ix]) py.show() # extract the betas beta = fit.extract('beta') ints_pa = pd.DataFrame(beta['beta'][:,:,0]) mpd_pa = pd.DataFrame(beta['beta'][:,:,1]) mtd_pa = pd.DataFrame(beta['beta'][:,:,2]) origin_pa = pd.DataFrame(beta['beta'][:,:,3]) mpdOrigin_pa = pd.DataFrame(beta['beta'][:,:,4]) mtdOrigin_pa = pd.DataFrame(beta['beta'][:,:,5]) # extract the kappas kappa = fit.extract('kappa') ints_rc = pd.DataFrame(kappa['kappa'][:,:,0]) mpd_rc = pd.DataFrame(kappa['kappa'][:,:,1]) mtd_rc = pd.DataFrame(kappa['kappa'][:,:,2]) origin_rc = pd.DataFrame(kappa['kappa'][:,:,3]) mpdOrigin_rc = pd.DataFrame(kappa['kappa'][:,:,4]) mtdOrigin_rc = pd.DataFrame(kappa['kappa'][:,:,5]) # check the plots DensPlot(ints_pa); TracePlot(ints_pa) DensPlot(mpd_pa); TracePlot(mpd_pa) DensPlot(mtd_pa); TracePlot(mtd_pa) DensPlot(origin_pa); TracePlot(origin_pa) DensPlot(mpdOrigin_pa); TracePlot(mpdOrigin_pa) DensPlot(mtdOrigin_pa); TracePlot(mtdOrigin_pa) DensPlot(ints_rc); TracePlot(ints_rc) DensPlot(mpd_rc); TracePlot(mpd_rc) DensPlot(mtd_rc); TracePlot(mtd_rc) DensPlot(origin_rc); TracePlot(origin_rc) DensPlot(mpdOrigin_rc); TracePlot(mpdOrigin_rc) DensPlot(mtdOrigin_rc); TracePlot(mtdOrigin_rc) ### PRESENCE/ABSENCE STAND LEVEL #### # first, make a color list cols = ['white', 'grey', 'black'] fig = py.figure(figsize = (4*5, 5.5)) # MPD-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mpd_pa.quantile(0.025)[x] > 0 or mpd_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mpd_pa.quantile(0.1)[x] > 0 or mpd_pa.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax1 = fig.add_subplot(151) ax1.hlines( envPlot['age'], mpd_pa.quantile(0.025), mpd_pa.quantile(0.975)) ax1.hlines( envPlot['age'], mpd_pa.quantile(0.1), mpd_pa.quantile(0.9), linewidth = 3) ax1.scatter( mpd_pa.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax1.axvline(0, linestyle = 'dashed', color = 'k') ax1.set_ylim([0, 200]) ax1.set_ylabel('Forest Age (Years)') ax1.set_title('MPD') # MTD-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mtd_pa.quantile(0.025)[x] > 0 or mtd_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtd_pa.quantile(0.1)[x] > 0 or mtd_pa.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax2 = fig.add_subplot(152, sharey = ax1) ax2.hlines( envPlot['age'], mtd_pa.quantile(0.025), mtd_pa.quantile(0.975)) ax2.hlines( envPlot['age'], mtd_pa.quantile(0.1), mtd_pa.quantile(0.9), linewidth = 3) ax2.scatter( mtd_pa.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax2.axvline(0, linestyle = 'dashed', color = 'k') ax2.set_title('MTD') ax2.set_ylim([0, 200]) py.setp(ax2.get_yticklabels(), visible = False) # ORIGIN-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(origin_pa.quantile(0.025)[x] > 0 or origin_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(origin_pa.quantile(0.1)[x] > 0 or origin_pa.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax3 = fig.add_subplot(153, sharey = ax1) ax3.hlines( envPlot['age'], origin_pa.quantile(0.025), origin_pa.quantile(0.975)) ax3.hlines( envPlot['age'], origin_pa.quantile(0.1), origin_pa.quantile(0.9), linewidth = 3) ax3.scatter( origin_pa.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax3.axvline(0, linestyle = 'dashed', color = 'k') ax3.set_title('Introduced') ax3.set_ylim([0, 200]) py.setp(ax3.get_yticklabels(), visible = False) # MPD:ORIGIN-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mpdOrigin_pa.quantile(0.025)[x] > 0 or mpdOrigin_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mpdOrigin_pa.quantile(0.1)[x] > 0 or mpdOrigin_pa.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax4 = fig.add_subplot(154, sharey = ax1) ax4.hlines( envPlot['age'], mpdOrigin_pa.quantile(0.025), mpdOrigin_pa.quantile(0.975)) ax4.hlines( envPlot['age'], mpdOrigin_pa.quantile(0.1), mpdOrigin_pa.quantile(0.9), linewidth = 3) ax4.scatter( mpdOrigin_pa.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax4.axvline(0, linestyle = 'dashed', color = 'k') ax4.set_title('MPD:Introduced') ax4.set_ylim([0, 200]) ax4.set_xticks([-2, -1, 0, 1, 2]) py.setp(ax4.get_yticklabels(), visible = False) # MTD:ORIGIN-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mtdOrigin_pa.quantile(0.025)[x] > 0 or mtdOrigin_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtdOrigin_pa.quantile(0.1)[x] > 0 or mtdOrigin_pa.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax5 = fig.add_subplot(155, sharey = ax1) ax5.hlines( envPlot['age'], mtdOrigin_pa.quantile(0.025), mtdOrigin_pa.quantile(0.975)) ax5.hlines( envPlot['age'], mtdOrigin_pa.quantile(0.1), mtdOrigin_pa.quantile(0.9), linewidth = 3) ax5.scatter( mtdOrigin_pa.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax5.axvline(0, linestyle = 'dashed', color = 'k') ax5.set_title('MTD:Introduced') ax5.set_ylim([0, 200]) ax5.set_ylabel('Presence/Absence') ax5.yaxis.set_label_position("right") py.setp(ax5.get_yticklabels(), visible = False) py.gca().invert_yaxis() py.savefig('STAN_ModelOut_PA_wtMPD.pdf', bbox_inches = 'tight') py.show() ## ENVIRONMENTAL PARAMETERS PLOT - PRESENCE ABSENCE f = py.figure(figsize = (4*5, 5.5)) # MPD ax1 = f.add_subplot(151) envMPD_PA = pd.DataFrame(fit.extract(['mpd_mu', 'mpd_light', 'mpd_vwc', 'mpd_litter' ,'mpd_age'])) sig95 = [int(envMPD_PA.quantile(0.025)[x] > 0 or envMPD_PA.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMPD_PA.quantile(0.1)[x] > 0 or envMPD_PA.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax1.scatter(envMPD_PA.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax1.hlines(range(5), envMPD_PA.quantile(0.025), envMPD_PA.quantile(0.975)) ax1.hlines(range(5), envMPD_PA.quantile(0.1), envMPD_PA.quantile(0.9), linewidth = 3) ax1.axvline(0, linestyle = 'dashed', color = 'k') ax1.set_ylim([-1, 5]) ax1.set_ylabel('Environmental Parameter') ax1.set_title('MPD') ax1.set_yticklabels(['Overall Effect', 'Light', 'VWC', 'Litter', 'Forest Age']) ax1.set_yticks(range(5)) ax1.set_xticks([-.4, -.2, 0, .2, .4]) # MTD ax2 = f.add_subplot(152, sharey = ax1) envMTD_PA = pd.DataFrame(fit.extract(['mtd_mu', 'mtd_light', 'mtd_vwc', 'mtd_litter', 'mtd_age'])) sig95 = [int(envMTD_PA.quantile(0.025)[x] > 0 or envMTD_PA.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMTD_PA.quantile(0.1)[x] > 0 or envMTD_PA.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax2.scatter(envMTD_PA.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax2.hlines(range(5), envMTD_PA.quantile(0.025), envMTD_PA.quantile(0.975)) ax2.hlines(range(5), envMTD_PA.quantile(0.1), envMTD_PA.quantile(0.9), linewidth = 3) ax2.axvline(0, linestyle = 'dashed', color = 'k') ax2.set_ylim([-1, 5]) ax2.set_title('MTD') ax2.set_yticks(range(5)) py.setp(ax2.get_yticklabels(), visible = False) # Origin ax3 = f.add_subplot(153, sharey = ax1) envOrigin_PA = pd.DataFrame(fit.extract(['origin_mu', 'origin_light', 'origin_vwc', 'origin_litter', 'origin_age'])) sig95 = [int(envOrigin_PA.quantile(0.025)[x] > 0 or envOrigin_PA.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envOrigin_PA.quantile(0.1)[x] > 0 or envOrigin_PA.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax3.scatter(envOrigin_PA.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax3.hlines(range(5), envOrigin_PA.quantile(0.025), envOrigin_PA.quantile(0.975)) ax3.hlines(range(5), envOrigin_PA.quantile(0.1), envOrigin_PA.quantile(0.9), linewidth = 3) ax3.axvline(0, linestyle = 'dashed', color = 'k') ax3.set_ylim([-1, 5]) ax3.set_title('Introduced') ax3.set_yticks(range(5)) ax3.set_xticks([-.8, -.4, 0, .4, .8]) py.setp(ax3.get_yticklabels(), visible = False) # MPD:Origin ax4 = f.add_subplot(154, sharey = ax1) envMPDOrigin_PA = pd.DataFrame(fit.extract(['mpdOrigin_mu', 'mpdOrigin_light', 'mpdOrigin_vwc', 'mpdOrigin_litter', 'mpdOrigin_age'])) sig95 = [int(envMPDOrigin_PA.quantile(0.025)[x] > 0 or envMPDOrigin_PA.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMPDOrigin_PA.quantile(0.1)[x] > 0 or envMPDOrigin_PA.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax4.scatter(envMPDOrigin_PA.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax4.hlines(range(5), envMPDOrigin_PA.quantile(0.025), envMPDOrigin_PA.quantile(0.975)) ax4.hlines(range(5), envMPDOrigin_PA.quantile(0.1), envMPDOrigin_PA.quantile(0.9), linewidth = 3) ax4.axvline(0, linestyle = 'dashed', color = 'k') ax4.set_ylim([-1, 5]) ax4.set_title('MPD:Introduced') ax4.set_yticks(range(5)) ax4.set_xticks([-0.4, 0, 0.4, 0.8]) py.setp(ax4.get_yticklabels(), visible = False) # MTD:Origin ax5 = f.add_subplot(155, sharey = ax1) envMTDOrigin_PA = pd.DataFrame(fit.extract(['mtdOrigin_mu', 'mtdOrigin_light', 'mtdOrigin_vwc', 'mtdOrigin_litter', 'mtdOrigin_age'])) sig95 = [int(envMTDOrigin_PA.quantile(0.025)[x] > 0 or envMTDOrigin_PA.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMTDOrigin_PA.quantile(0.1)[x] > 0 or envMTDOrigin_PA.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax5.scatter(envMTDOrigin_PA.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax5.hlines(range(5), envMTDOrigin_PA.quantile(0.025), envMTDOrigin_PA.quantile(0.975)) ax5.hlines(range(5), envMTDOrigin_PA.quantile(0.1), envMTDOrigin_PA.quantile(0.9), linewidth = 3) ax5.axvline(0, linestyle = 'dashed', color = 'k') ax5.set_ylim([-1, 5]) ax5.set_title('MPD:Introduced') ax5.set_yticks(range(5)) ax5.set_xticks([-.8, -.4, 0, .4, .8]) ax5.set_ylabel('Presence/Absence') ax5.yaxis.set_label_position("right") py.setp(ax5.get_yticklabels(), visible = False) py.gca().invert_yaxis() py.savefig('envQuants_PA_STAN_wtMPD.pdf', bbox_inches = 'tight') py.show() ### RELATIVE COVER STAND LEVEL COEFFICIENTS ### fig = py.figure(figsize = (4*5, 5.5)) # MPD-PA # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mpd_rc.quantile(0.025)[x] > 0 or mpd_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mpd_rc.quantile(0.1)[x] > 0 or mpd_rc.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax1 = fig.add_subplot(151) ax1.hlines( envPlot['age'], mpd_rc.quantile(0.025), mpd_rc.quantile(0.975)) ax1.hlines( envPlot['age'], mpd_rc.quantile(0.1), mpd_rc.quantile(0.9), linewidth = 3) ax1.scatter( mpd_rc.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax1.axvline(0, linestyle = 'dashed', color = 'k') ax1.set_ylim([0, 200]) ax1.set_xticks([-2,-1, 0, 1, 2]) ax1.set_ylabel('Forest Age (Years)') ax1.set_title('MPD') # MTD-RC # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mtd_rc.quantile(0.025)[x] > 0 or mtd_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtd_rc.quantile(0.1)[x] > 0 or mtd_rc.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax2 = fig.add_subplot(152, sharey = ax1) ax2.hlines( envPlot['age'], mtd_rc.quantile(0.025), mtd_rc.quantile(0.975)) ax2.hlines( envPlot['age'], mtd_rc.quantile(0.1), mtd_rc.quantile(0.9), linewidth = 3) ax2.scatter( mtd_rc.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax2.axvline(0, linestyle = 'dashed', color = 'k') ax2.set_title('MTD') ax2.set_ylim([0, 200]) ax2.set_xticks([-2, -1, 0, 1, 2]) py.setp(ax2.get_yticklabels(), visible = False) # ORIGIN-RC # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(origin_rc.quantile(0.025)[x] > 0 or origin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(origin_rc.quantile(0.1)[x] > 0 or origin_rc.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax3 = fig.add_subplot(153, sharey = ax1) ax3.hlines( envPlot['age'], origin_rc.quantile(0.025), origin_rc.quantile(0.975)) ax3.hlines( envPlot['age'], origin_rc.quantile(0.1), origin_rc.quantile(0.9), linewidth = 3) ax3.scatter( origin_rc.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax3.axvline(0, linestyle = 'dashed', color = 'k') ax3.set_title('Introduced') ax3.set_ylim([0, 200]) py.setp(ax3.get_yticklabels(), visible = False) # MPD:ORIGIN-RC # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mpdOrigin_rc.quantile(0.025)[x] > 0 or mpdOrigin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mpdOrigin_rc.quantile(0.1)[x] > 0 or mpdOrigin_rc.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax4 = fig.add_subplot(154, sharey = ax1) ax4.hlines( envPlot['age'], mpdOrigin_rc.quantile(0.025), mpdOrigin_rc.quantile(0.975)) ax4.hlines( envPlot['age'], mpdOrigin_rc.quantile(0.1), mpdOrigin_rc.quantile(0.9), linewidth = 3) ax4.scatter( mpdOrigin_rc.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax4.axvline(0, linestyle = 'dashed', color = 'k') ax4.set_title('MPD:Introduced') ax4.set_ylim([0, 200]) ax4.set_xticks([-2, -1, 0, 1, 2]) py.setp(ax4.get_yticklabels(), visible = False) # MTD:ORIGIN-RC # first, determine if the parameter is significant at the 95% quantile level sig95 = [int(mtdOrigin_rc.quantile(0.025)[x] > 0 or mtdOrigin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtdOrigin_rc.quantile(0.1)[x] > 0 or mtdOrigin_rc.quantile(0.9)[x] < 0) for x in range(25)] # add the lists together so that 0 indicates no significance, a 1 indicates 80%, and a 2 indicates 95% sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax5 = fig.add_subplot(155, sharey = ax1) ax5.hlines( envPlot['age'], mtdOrigin_rc.quantile(0.025), mtdOrigin_rc.quantile(0.975)) ax5.hlines( envPlot['age'], mtdOrigin_rc.quantile(0.1), mtdOrigin_rc.quantile(0.9), linewidth = 3) ax5.scatter( mtdOrigin_rc.median(), envPlot['age'], s = 100, color = sigCols, edgecolor = 'k', zorder = 10) ax5.axvline(0, linestyle = 'dashed', color = 'k') ax5.set_title('MTD:Introduced') ax5.set_ylim([0, 200]) ax5.set_ylabel('Relative Cover') ax5.yaxis.set_label_position("right") py.setp(ax5.get_yticklabels(), visible = False) py.gca().invert_yaxis() py.savefig('STAN_ModelOut_RC_wtMPD.pdf', bbox_inches = 'tight') py.show() ## ENVIRONMENTAL PARAMETERS PLOT - RELATIVE COVER f = py.figure(figsize = (4*5, 5.5)) # MPD ax1 = f.add_subplot(151) envMPD_RC = pd.DataFrame(fit.extract(['mpd_mu_RC', 'mpd_light_RC', 'mpd_vwc_RC', 'mpd_litter_RC', 'mpd_age_RC'])) sig95 = [int(envMPD_RC.quantile(0.025)[x] > 0 or envMPD_RC.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMPD_RC.quantile(0.1)[x] > 0 or envMPD_RC.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax1.scatter(envMPD_RC.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax1.hlines(range(5), envMPD_RC.quantile(0.025), envMPD_RC.quantile(0.975)) ax1.hlines(range(5), envMPD_RC.quantile(0.1), envMPD_RC.quantile(0.9), linewidth = 3) ax1.axvline(0, linestyle = 'dashed', color = 'k') ax1.set_ylim([-1, 5]) ax1.set_ylabel('Environmental Parameter') ax1.set_title('MPD') ax1.set_yticklabels(['Overall Effect', 'Light', 'VWC', 'Litter', 'Forest Age']) ax1.set_yticks(range(5)) ax1.set_xticks([-.8, -.4, 0, .4, .8]) # MTD ax2 = f.add_subplot(152, sharey = ax1) envMTD_RC = pd.DataFrame(fit.extract(['mtd_mu_RC', 'mtd_light_RC', 'mtd_vwc_RC', 'mtd_litter_RC', 'mtd_age_RC'])) sig95 = [int(envMTD_RC.quantile(0.025)[x] > 0 or envMTD_RC.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMTD_RC.quantile(0.1)[x] > 0 or envMTD_RC.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax2.scatter(envMTD_RC.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax2.hlines(range(5), envMTD_RC.quantile(0.025), envMTD_RC.quantile(0.975)) ax2.hlines(range(5), envMTD_RC.quantile(0.1), envMTD_RC.quantile(0.9), linewidth = 3) ax2.axvline(0, linestyle = 'dashed', color = 'k') ax2.set_ylim([-1, 5]) ax2.set_title('MTD') ax2.set_yticks(range(5)) ax2.set_xticks([-.8, -.4, 0, .4, .8]) py.setp(ax2.get_yticklabels(), visible = False) # Origin ax3 = f.add_subplot(153, sharey = ax1) envOrigin_RC = pd.DataFrame(fit.extract(['origin_mu_RC', 'origin_light_RC', 'origin_vwc_RC', 'origin_litter_RC', 'origin_age_RC'])) sig95 = [int(envOrigin_RC.quantile(0.025)[x] > 0 or envOrigin_RC.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envOrigin_RC.quantile(0.1)[x] > 0 or envOrigin_RC.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax3.scatter(envOrigin_RC.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax3.hlines(range(5), envOrigin_RC.quantile(0.025), envOrigin_RC.quantile(0.975)) ax3.hlines(range(5), envOrigin_RC.quantile(0.1), envOrigin_RC.quantile(0.9), linewidth = 3) ax3.axvline(0, linestyle = 'dashed', color = 'k') ax3.set_ylim([-1, 5]) ax3.set_title('Introduced') ax3.set_yticks(range(5)) ax3.set_xticks([-.8, -.4, 0, .4, .8]) py.setp(ax3.get_yticklabels(), visible = False) # MPD:Origin ax4 = f.add_subplot(154, sharey = ax1) envMPDOrigin_RC = pd.DataFrame(fit.extract(['mpdOrigin_mu_RC', 'mpdOrigin_light_RC', 'mpdOrigin_vwc_RC', 'mpdOrigin_litter_RC', 'mpdOrigin_age_RC'])) sig95 = [int(envMPDOrigin_RC.quantile(0.025)[x] > 0 or envMPDOrigin_RC.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMPDOrigin_RC.quantile(0.1)[x] > 0 or envMPDOrigin_RC.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax4.scatter(envMPDOrigin_RC.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax4.hlines(range(5), envMPDOrigin_RC.quantile(0.025), envMPDOrigin_RC.quantile(0.975)) ax4.hlines(range(5), envMPDOrigin_RC.quantile(0.1), envMPDOrigin_RC.quantile(0.9), linewidth = 3) ax4.axvline(0, linestyle = 'dashed', color = 'k') ax4.set_ylim([-1, 5]) ax4.set_title('MPD:Introduced') ax4.set_yticks(range(5)) ax4.set_xticks([-.8, -.4, 0, .4, .8]) py.setp(ax4.get_yticklabels(), visible = False) # MTD:Origin ax5 = f.add_subplot(155, sharey = ax1) envMTDOrigin_RC = pd.DataFrame(fit.extract(['mtdOrigin_mu_RC', 'mtdOrigin_light_RC', 'mtdOrigin_vwc_RC', 'mtdOrigin_litter_RC', 'mtdOrigin_age_RC'])) sig95 = [int(envMTDOrigin_RC.quantile(0.025)[x] > 0 or envMTDOrigin_RC.quantile(0.975)[x] < 0) for x in range(5)] sig80 = [int(envMTDOrigin_RC.quantile(0.1)[x] > 0 or envMTDOrigin_RC.quantile(0.9)[x] < 0) for x in range(5)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] ax5.scatter(envMTDOrigin_RC.median(), range(5), color = sigCols, s = 100, zorder = 10, edgecolor = 'k') ax5.hlines(range(5), envMTDOrigin_RC.quantile(0.025), envMTDOrigin_RC.quantile(0.975)) ax5.hlines(range(5), envMTDOrigin_RC.quantile(0.1), envMTDOrigin_RC.quantile(0.9), linewidth = 3) ax5.axvline(0, linestyle = 'dashed', color = 'k') ax5.set_ylim([-1, 5]) ax5.set_title('MPD:Introduced') ax5.set_yticks(range(5)) ax5.set_ylabel('Relative Cover') ax5.yaxis.set_label_position("right") py.setp(ax5.get_yticklabels(), visible = False) py.gca().invert_yaxis() py.savefig('envQuants_RC_STAN_wtMPD.pdf', bbox_inches = 'tight') py.show() ### NEXT, INDIVIDUAL FOREST PREDICTIONS # put the forest in order of age plotAge = envPlot[['plot', 'age']].sort('age') plotAge['title'] = plotAge.apply(lambda x: '%i (%s)' % (x['age'], x['plot']), axis = 1) # PRESENCE-ABSENCE : MPD # make the figure f = py.figure(figsize = (20, 20)) # for each forest.. for i in range(25): ix = plotAge['age'].index[i] # get the index value associate with the first aged plot plotNames = plotAge['plot'][ix] ax = f.add_subplot(5, 5, i+1) # make the axis mpdPred = np.linspace(min(fullData[fullData['plot'] == plotNames]['MPD_scale']) , max(fullData[fullData['plot'] == plotNames]['MPD_scale']) ) # get MPD scale for the plot nPred = ints_pa.median()[ix] + mpd_pa.median()[ix]*mpdPred # get the prediction for native species nPred = np.exp(nPred) / (1 + np.exp(nPred)) # NATIVE SPECIES FIRST if mpd_pa.quantile(0.025)[ix] > 0 or mpd_pa.quantile(0.975)[ix] < 0: ax.plot(mpdPred, nPred, 'k-', linewidth = 2) elif mpd_pa.quantile(0.1)[ix] > 0 or mpd_pa.quantile(0.9)[ix] < 0: ax.plot(mpdPred, nPred, 'k', linestyle = 'dashed', linewidth = 2) else: ax.plot(mpdPred, nPred, 'k', linestyle = 'dotted', linewidth = 2) # IF THE INTERACTION IS EVEN MODERATLEY SIGNIFICANT, CALCULATE THE NEW SLOPE if mpdOrigin_pa.quantile(0.1)[ix] > 0 or mpdOrigin_pa.quantile(0.9)[ix] < 0: intSlope = mpd_pa[ix] + mpdOrigin_pa[ix] # IF THE NEW SLOPE IS SIGNIFICANT, PLOT if intSlope.quantile(0.025) > 0 or intSlope.quantile(0.975) < 0: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mpdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mpdPred, iPred, 'r-', linewidth = 2) elif intSlope.quantile(0.1) > 0 or intSlope.quantile(0.9) < 0: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mpdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mpdPred, iPred, 'r--', linewidth = 2) else: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mpdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mpdPred, iPred, 'r', linestyle = 'dotted', linewidth = 2) ax.set_ylim([0,1]) ax.set_title(plotAge['title'][ix]) ax.set_xticks(range(-4, 3)) f.text(0.5, 0.075, 'Standardized MPD', va = 'center', ha = 'center', fontsize = 14) f.text(0.1, 0.5, 'Probability of Occurrence', va = 'center', ha = 'center', fontsize = 14, rotation = 90) py.savefig('MPD_PA_stands_wtMPD.pdf', bbox_inches = 'tight') py.show() # PRESENCE-ABSENCE : MTD # make the figure f = py.figure(figsize = (20, 20)) # for each forest.. for i in range(25): ix = plotAge['age'].index[i] # get the index value associate with the first aged plot plotNames = plotAge['plot'][ix] ax = f.add_subplot(5, 5, i+1) # make the axis mtdPred = np.linspace(min(fullData[fullData['plot'] == plotNames]['MTD_scale']) , max(fullData[fullData['plot'] == plotNames]['MTD_scale']) ) # get MTD scale for the plot nPred = ints_pa.median()[ix] + mtd_pa.median()[ix]*mtdPred # get the prediction for native species nPred = np.exp(nPred) / (1 + np.exp(nPred)) # NATIVE SPECIES FIRST if mtd_pa.quantile(0.025)[ix] > 0 or mtd_pa.quantile(0.975)[ix] < 0: ax.plot(mtdPred, nPred, 'k-', linewidth = 2) elif mtd_pa.quantile(0.1)[ix] > 0 or mtd_pa.quantile(0.9)[ix] < 0: ax.plot(mtdPred, nPred, 'k', linestyle = 'dashed', linewidth = 2) else: ax.plot(mtdPred, nPred, 'k', linestyle = 'dotted', linewidth = 2) # IF THE INTERACTION IS EVEN MODERATLEY SIGNIFICANT, CALCULATE THE NEW SLOPE if mtdOrigin_pa.quantile(0.1)[ix] > 0 or mtdOrigin_pa.quantile(0.9)[ix] < 0: intSlope = mtd_pa[ix] + mtdOrigin_pa[ix] # IF THE NEW SLOPE IS SIGNIFICANT, PLOT if intSlope.quantile(0.025) > 0 or intSlope.quantile(0.975) < 0: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mtdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mtdPred, iPred, 'r-', linewidth = 2) elif intSlope.quantile(0.1) > 0 or intSlope.quantile(0.9) < 0: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mtdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mtdPred, iPred, 'r--', linewidth = 2) else: iPred = ints_pa.median()[ix] + origin_pa.median()[ix] + intSlope.median()*mtdPred iPred = np.exp(iPred) / (1 + np.exp(iPred)) ax.plot(mtdPred, iPred, 'r', linestyle = 'dotted', linewidth = 2) ax.set_ylim([0,1]) ax.set_title(plotAge['title'][ix]) ax.set_xticks(range(-2, 5)) f.text(0.5, 0.075, 'Standardized MTD', va = 'center', ha = 'center', fontsize = 14) f.text(0.1, 0.5, 'Probability of Occurrence', va = 'center', ha = 'center', fontsize = 14, rotation = 90) py.savefig('MTD_PA_stands_wtMPD.pdf', bbox_inches = 'tight') py.show() # RELATIVE COVER : MPD # make the figure f = py.figure(figsize = (20, 20)) # for each forest.. for i in range(25): ix = plotAge['age'].index[i] # get the index value associate with the first aged plot plotNames = plotAge['plot'][ix] ax = f.add_subplot(5, 5, i+1) # make the axis mpdPred = np.linspace(min(fullData[fullData['plot'] == plotNames]['MPD_scale']) , max(fullData[fullData['plot'] == plotNames]['MPD_scale']) ) # get MPD scale for the plot nPred = ints_rc.median()[ix] + mpd_rc.median()[ix]*mpdPred # get the prediction for native species nPred = np.exp(nPred) # NATIVE SPECIES FIRST if mpd_rc.quantile(0.025)[ix] > 0 or mpd_rc.quantile(0.975)[ix] < 0: ax.plot(mpdPred, nPred, 'k-', linewidth = 2) elif mpd_rc.quantile(0.1)[ix] > 0 or mpd_rc.quantile(0.9)[ix] < 0: ax.plot(mpdPred, nPred, 'k', linestyle = 'dashed', linewidth = 2) else: ax.plot(mpdPred, nPred, 'k', linestyle = 'dotted', linewidth = 2) # IF THE INTERACTION IS EVEN MODERATLEY SIGNIFICANT, CALCULATE THE NEW SLOPE if mpdOrigin_rc.quantile(0.1)[ix] > 0 or mpdOrigin_rc.quantile(0.9)[ix] < 0: intSlope = mpd_rc[ix] + mpdOrigin_rc[ix] iPred = ints_rc.median()[ix] + origin_rc.median()[ix] + intSlope.median()*mtdPred iPred = np.exp(iPred) # IF THE NEW SLOPE IS SIGNIFICANT, PLOT if intSlope.quantile(0.025) > 0 or intSlope.quantile(0.975) < 0: ax.plot(mpdPred, iPred, 'r-', linewidth = 2) elif intSlope.quantile(0.1) > 0 or intSlope.quantile(0.9) < 0: ax.plot(mpdPred, iPred, 'r--', linewidth = 2) else: ax.plot(mpdPred, iPred, 'r', linestyle = 'dotted', linewidth = 2) ax.set_ylim([0,0.5]) ax.set_title(plotAge['title'][ix]) ax.set_xticks(range(-4, 3)) f.text(0.5, 0.075, 'Standardized MPD', va = 'center', ha = 'center', fontsize = 14) f.text(0.1, 0.5, 'Relative Cover', va = 'center', ha = 'center', fontsize = 14, rotation = 90) py.savefig('MPD_RC_stands_wtMPD.pdf', bbox_inches = 'tight') py.show() # RELATIVE COVER : MTD # make the figure f = py.figure(figsize = (20, 20)) # for each forest.. for i in range(25): ix = plotAge['age'].index[i] # get the index value associate with the first aged plot plotNames = plotAge['plot'][ix] ax = f.add_subplot(5, 5, i+1) # make the axis mtdPred = np.linspace(min(fullData[fullData['plot'] == plotNames]['MTD_scale']) , max(fullData[fullData['plot'] == plotNames]['MTD_scale']) ) # get MTD scale for the plot nPred = ints_rc.median()[ix] + mtd_rc.median()[ix]*mtdPred # get the prediction for native species nPred = np.exp(nPred) # NATIVE SPECIES FIRST if mtd_rc.quantile(0.025)[ix] > 0 or mtd_rc.quantile(0.975)[ix] < 0: ax.plot(mtdPred, nPred, 'k-', linewidth = 2) elif mtd_rc.quantile(0.1)[ix] > 0 or mtd_rc.quantile(0.9)[ix] < 0: ax.plot(mtdPred, nPred, 'k', linestyle = 'dashed', linewidth = 2) else: ax.plot(mtdPred, nPred, 'k', linestyle = 'dotted', linewidth = 2) # IF THE INTERACTION IS EVEN MODERATLEY SIGNIFICANT, CALCULATE THE NEW SLOPE if mtdOrigin_rc.quantile(0.1)[ix] > 0 or mtdOrigin_rc.quantile(0.9)[ix] < 0: intSlope = mtd_rc[ix] + mtdOrigin_rc[ix] iPred = ints_rc.median()[ix] + origin_rc.median()[ix] + intSlope.median()*mtdPred iPred = np.exp(iPred) # IF THE NEW SLOPE IS SIGNIFICANT, PLOT if intSlope.quantile(0.025) > 0 or intSlope.quantile(0.975) < 0: ax.plot(mtdPred, iPred, 'r-', linewidth = 2) elif intSlope.quantile(0.1) > 0 or intSlope.quantile(0.9) < 0: ax.plot(mtdPred, iPred, 'r--', linewidth = 2) else: ax.plot(mtdPred, iPred, 'r', linestyle = 'dotted', linewidth = 2) ax.set_ylim([0,0.1]) ax.set_title(plotAge['title'][ix]) ax.set_xticks(range(-2, 5)) f.text(0.5, 0.075, 'Standardized MTD', va = 'center', ha = 'center', fontsize = 14) f.text(0.1, 0.5, 'Relative Cover', va = 'center', ha = 'center', fontsize = 14, rotation = 90) py.savefig('MTD_RC_stands_wtMPD.pdf', bbox_inches = 'tight') py.show() # ALL MPD IN ONE PLOT mpdPred = np.linspace(min(fullData['MPD_scale']), max(fullData['MPD_scale']), num = 100) from matplotlib.lines import Line2D for i in range(25): nPred = np.exp(ints_rc.median()[i] + mpd_rc.median()[i]*mpdPred) if mpd_rc.quantile(0.025)[i] > 0 or mpd_rc.quantile(0.975)[i] < 0: py.plot(mpdPred, nPred, 'k', linestyle = 'solid') elif mpd_rc.quantile(0.1)[i] > 0 or mpd_rc.quantile(0.9)[i] < 0: py.plot(mpdPred, nPred, 'k', linestyle = 'dashed') else: py.plot(mpdPred, nPred, 'k', linestyle = 'dotted') l1 = Line2D([], [], linewidth=1, color = 'k', linestyle="dotted") l2 = Line2D([], [], linewidth=1, color = 'k', linestyle="dashed") l3 = Line2D([], [], linewidth=1, color = 'k', linestyle="solid") py.legend([l1, l2, l3], ["Not Significant", "80% CI", "95% CI"], loc = 1) py.xlabel('Standardized MPD') py.ylabel('Relative Cover of Native Species') py.savefig('MPD_SingleFig_wtMPD.pdf', bbox_inches = 'tight') py.show() ############################# ############################# ### NEXT, GROUP LEVEL EFFECTS from matplotlib.patches import Polygon ## LITTER - INTRODUCED PA origin_int = fit.extract('origin_mu') origin_litter = fit.extract('origin_litter') sig95 = [int(origin_pa.quantile(0.025)[x] > 0 or origin_pa.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(origin_pa.quantile(0.1)[x] > 0 or origin_pa.quantile(0.9)[x] < 0) for x in range(25)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] litterPred = np.array([min(envPlot['litter_scale']), max(envPlot['litter_scale'])]) yhat = np.median(origin_int['origin_mu']) + np.median(origin_litter['origin_litter'])*litterPred f = py.figure() ax = f.add_subplot(111) ax.scatter(envPlot['litter_scale'], origin_pa.median(), c = sigCols, s = 150, zorder = 10) ax.plot(litterPred, yhat, 'k-') ax.add_patch(Polygon([[-1.5, -2], [2.5, -2], [2.5, 0], [-1.5, 0]], closed = True, fill = True, fc = 'lightgrey', linestyle = 'dashed')) ax.text(0.5, -1.75, 'Introduced Species\nLess Likely to Occur', ha = 'center', va = 'center') ax.text(1.5, 0.25, 'Introduced Species\nMore Likely to Occur', ha = 'center', va = 'center') py.ylim([-2, 1]) py.xlim([-1.5, 2.5]) py.xlabel('Standardized Litter Depth') py.ylabel('Introduced Presence/Absence Coefficient') py.savefig('Origin_Litter.pdf', bbox_inches = 'tight') py.show() ## WATER - INTRODUCED RC oINT = fit.extract('origin_mu_RC') oVWC = fit.extract('origin_vwc_RC') sig95 = [int(origin_rc.quantile(0.025)[x] > 0 or origin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(origin_rc.quantile(0.1)[x] > 0 or origin_rc.quantile(0.9)[x] < 0) for x in range(25)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] vwcPred = np.array([min(envPlot['water_scale']), max(envPlot['water_scale'])]) yhat = np.median(oINT['origin_mu_RC']) + np.median(oVWC['origin_vwc_RC'])*vwcPred f = py.figure() ax = f.add_subplot(111) ax.scatter(envPlot['water_scale'], origin_rc.median(), c = sigCols, s = 150, zorder = 10) ax.plot(vwcPred, yhat, 'k-') ax.add_patch(Polygon([[-2.2, -2], [2.6, -2], [2.6, 0], [-2.2, 0]], closed = True, fill = True, fc = 'lightgrey', linestyle = 'dashed')) ax.text(0.5, -1.25, 'Introduced Species\nLess Abundant Than Natives', ha = 'center', va = 'center') ax.text(1.4, 1.3, 'Introduced Species\nMore Abundant Than Natives', ha = 'center', va = 'center') py.ylim([-1.5, 1.5]) py.xlim([-2.2, 2.6]) py.xlabel('Standardized VWC') py.ylabel('Introduced Relative Cover Coefficient') py.savefig('Origin_Water.pdf', bbox_inches = 'tight') py.show() np.mean(oVWC['origin_vwc_RC'] > 0) ## FOREST AGE - INTRODUCED RC oINT = fit.extract('origin_mu_RC') oAge = fit.extract('origin_age_RC') sig95 = [int(origin_rc.quantile(0.025)[x] > 0 or origin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(origin_rc.quantile(0.1)[x] > 0 or origin_rc.quantile(0.9)[x] < 0) for x in range(25)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] agePred = np.array([min(envPlot['age_scale']), max(envPlot['age_scale'])]) yhat = np.median(oINT['origin_mu_RC']) + np.median(oAge['origin_age_RC'])*agePred f = py.figure() ax = f.add_subplot(111) ax.scatter(envPlot['age_scale'], origin_rc.median(), c = sigCols, s = 150, zorder = 10) ax.plot(vwcPred, yhat, 'k-') ax.add_patch(Polygon([[-2.2, -2], [2.6, -2], [2.6, 0], [-2.2, 0]], closed = True, fill = True, fc = 'lightgrey', linestyle = 'dashed')) ax.text(0.5, -1.25, 'Introduced Species\nLess Abundant Than Natives', ha = 'center', va = 'center') ax.text(1.4, 1.3, 'Introduced Species\nMore Abundant Than Natives', ha = 'center', va = 'center') py.ylim([-1.5, 1.5]) py.xlim([-2.2, 2.6]) py.xlabel('Standardized Forest Age') py.ylabel('Introduced Relative Cover Coefficient') py.savefig('Origin_Age.pdf', bbox_inches = 'tight') py.show() np.mean(oAge['origin_age_RC'] < 0) ## LIGHT - MTD:ORIGIN RC mtdOrigin_INT = fit.extract('mtdOrigin_mu_RC') mtdOrigin_Light = fit.extract('mtdOrigin_light_RC') sig95 = [int(mtdOrigin_rc.quantile(0.025)[x] > 0 or mtdOrigin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtdOrigin_rc.quantile(0.1)[x] > 0 or mtdOrigin_rc.quantile(0.9)[x] < 0) for x in range(25)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] lightPred = np.array([min(envPlot['light_scale']), max(envPlot['light_scale'])]) yhat = np.median(mtdOrigin_INT['mtdOrigin_mu_RC']) + np.median(mtdOrigin_Light['mtdOrigin_light_RC'])*lightPred f = py.figure() ax = f.add_subplot(111) ax.scatter(envPlot['light_scale'], mtdOrigin_rc.median(), c = sigCols, s = 150, zorder = 10) ax.plot(lightPred, yhat, 'k-') ax.add_patch(Polygon([[-2.2, -2], [2.6, -2], [2.6, 0], [-2.2, 0]], closed = True, fill = True, fc = 'lightgrey', linestyle = 'dashed')) ax.text(0.5, -0.75, 'Stronger Trait Filtering\nThan Natives', ha = 'center', va = 'center') ax.text(-1, 1.27, 'Weaker Trait Filtering\nThan Natives', ha = 'center', va = 'center') py.xlabel('Standardized Light Transmittance') py.ylabel('MTD:Introduced Relative Cover Coefficient') py.ylim([-1, 1.5]) py.xlim([-2.2, 2.6]) py.savefig('mtdOrigin_Light.pdf', bbox_inches = 'tight') py.show() ## AGE - MTD:ORIGIN RC mtdOrigin_INT = fit.extract('mtdOrigin_mu_RC') mtdOrigin_Age = fit.extract('mtdOrigin_age_RC') sig95 = [int(mtdOrigin_rc.quantile(0.025)[x] > 0 or mtdOrigin_rc.quantile(0.975)[x] < 0) for x in range(25)] sig80 = [int(mtdOrigin_rc.quantile(0.1)[x] > 0 or mtdOrigin_rc.quantile(0.9)[x] < 0) for x in range(25)] sigLevel = [sum(x) for x in zip(sig95, sig80)] sigCols = [cols[x] for x in sigLevel] agePred = np.array([min(envPlot['age_scale']), max(envPlot['age_scale'])]) yhat = np.median(mtdOrigin_INT['mtdOrigin_mu_RC']) + np.median(mtdOrigin_Age['mtdOrigin_age_RC'])*lightPred f = py.figure() ax = f.add_subplot(111) ax.scatter(envPlot['age_scale'], mtdOrigin_rc.median(), c = sigCols, s = 150, zorder = 10) ax.plot(agePred, yhat, 'k-') ax.add_patch(Polygon([[-2.2, -2], [2.6, -2], [2.6, 0], [-2.2, 0]], closed = True, fill = True, fc = 'lightgrey', linestyle = 'dashed')) ax.text(0.3, -1.5, 'Stronger Trait Filtering\nThan Natives', ha = 'center', va = 'center') ax.text(-0.5, 2, 'Weaker Trait Filtering\nThan Natives', ha = 'center', va = 'center') py.xlabel('Standardized Forest Age') py.ylabel('MTD:Introduced Relative Cover Coefficient') #py.ylim([-1, 1.5]) py.xlim([-2.2, 2.6]) py.savefig('mtdOrigin_Age.pdf', bbox_inches = 'tight') py.show()