"""Biological introductions as a large-sample base rate for 'arrives, establishes, spreads'. Three public datasets: 1. Moulton & Cropper 2019 (PeerJ 7:e7637) data file S1: 832 bird introduction events with number released and fate, originally compiled by Sol et al. 2012 (Science 337, 580). 2. GAVIA (Dyer, Redding & Blackburn 2017, Sci. Data 4:170041): 27,723 distribution records of alien birds. 3. Cassey et al. 2018 (PLoS Biol. 16:e2005987) S1 Data: effect sizes of the propagule-pressure meta-analysis. """ import io import os import zipfile import numpy as np import pandas as pd from common import DATA, SEED, dump_json, fetch, jeffreys, logistic_fit def propagule(): zpath = fetch("moulton_zip") with zipfile.ZipFile(zpath) as z: name = [n for n in z.namelist() if n.endswith(".csv")][0] d = pd.read_csv(io.BytesIO(z.read(name)), encoding="latin1") d = d[d["fate"].notna()].copy() d = d.rename(columns={d.columns[0]: "species"}) d["Number"] = d["Number"].astype(float) bins = [(2, 10), (11, 30), (31, 100), (101, 300), (301, 1000), (1001, 1e12)] rows = [] for lo, hi in bins: m = (d.Number >= lo) & (d.Number <= hi) k, n = int(d.fate[m].sum()), int(m.sum()) ci = jeffreys(k, n) rows.append(dict(bin=f"{lo}-{int(hi) if hi < 1e11 else 'more'}", lo=lo, hi=min(hi, 1e7), n=n, established=k, fraction=k / n, ci_lo=ci[0], ci_hi=ci[1], geo_mid=float(np.exp(np.log(d.Number[m]).mean())))) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(DATA, "propagule_establishment_bins.csv"), index=False) x = np.log10(d.Number.values) y = d.fate.values.astype(float) b, cov, ll = logistic_fit(x[:, None], y) b0, _, ll0 = logistic_fit(np.zeros((len(y), 0)), y) m3 = d.Number <= 1000 # game-bird mass releases above 1000 behave differently b3, cov3, _ = logistic_fit(x[m3.values][:, None], y[m3.values]) rng = np.random.default_rng(SEED + 3) boot = [] for _ in range(1000): j = rng.integers(0, len(y), len(y)) bb, _, _ = logistic_fit(x[j][:, None], y[j]) boot.append(bb) boot = np.array(boot) pred = lambda N, c=b: float(1 / (1 + np.exp(-(c[0] + c[1] * np.log10(N))))) predci = lambda N: [float(np.percentile(1 / (1 + np.exp(-(boot[:, 0] + boot[:, 1] * np.log10(N)))), q)) for q in (2.5, 97.5)] # species-clustered check: one event per species drawn at random, repeated sl = [] for _ in range(300): s = d.groupby("species").sample(1, random_state=int(rng.integers(1e9))) bb, _, _ = logistic_fit(np.log10(s.Number.values)[:, None], s.fate.values.astype(float)) sl.append(bb[1]) # how concentrated are the events in a few species? (species-level pseudo-replication) cnt = d.species.value_counts() cum = cnt.cumsum() / cnt.sum() n_half = int((cum < 0.5).sum() + 1) out = dict( source="Moulton & Cropper 2019 data file S1 (Sol et al. 2012 events)", n_events=int(len(d)), n_species=int(d.species.nunique()), overall_established=int(y.sum()), overall_fraction=float(y.mean()), overall_jeffreys95=jeffreys(int(y.sum()), len(y)), median_number_released=float(d.Number.median()), species_concentration=dict(species_supplying_half_of_events=n_half, share_of_top_23_species=float(cnt.iloc[:23].sum() / cnt.sum()), max_events_one_species=int(cnt.iloc[0])), literature_note="The authors of this data file concluded that propagule pressure does NOT consistently predict establishment " "(Moulton & Cropper 2019, PeerJ 7:e7637). Blackburn, Cassey, Lockwood & Duncan 2020 (PeerJ 8:e8766) re-analysed " "the same data and found a strong positive relation; the pooled slope fitted here agrees with the re-analysis.", bins=rows, check_against_published=dict(frac_2_10=float(d.fate[(d.Number >= 2) & (d.Number <= 10)].mean()), frac_11_100=float(d.fate[(d.Number >= 11) & (d.Number <= 100)].mean()), published_2_10=0.28, published_11_100=0.34), logistic_all=dict(intercept=float(b[0]), slope_per_dex=float(b[1]), se_slope=float(np.sqrt(cov[1, 1])), lr_chi2=float(2 * (ll - ll0)), boot95_slope=[float(np.percentile(boot[:, 1], 2.5)), float(np.percentile(boot[:, 1], 97.5))], one_event_per_species_slope_median=float(np.median(sl)), one_event_per_species_slope_p2p5_p97p5=[float(np.percentile(sl, 2.5)), float(np.percentile(sl, 97.5))]), logistic_upto1000=dict(intercept=float(b3[0]), slope_per_dex=float(b3[1]), se_slope=float(np.sqrt(cov3[1, 1])), n=int(m3.sum())), predicted_establishment={str(N): dict(p=pred(N), boot95=predci(N), p_upto1000_fit=pred(N, b3)) for N in (2, 10, 100, 1000)}, ) return out, tab, (b, b3) def gavia(): d = pd.read_csv(fetch("gavia"), low_memory=False, encoding="latin1") d["fail"] = d.StatusCat.isin(["Unsuccessful", "DiedOut", "Extirpated"]) d["est"] = d.StatusCat.eq("Established") d["known"] = d.fail | d.est d["area"] = d.AreaName1.fillna("") mode = lambda s: s.mode().iat[0] if s.notna().any() else np.nan g = d.groupby(["Binomial", "CountryName"]).agg(est=("est", "any"), known=("known", "any"), land=("LandType", mode)) k = g[g.known] by_land = {} for land, sub in k.groupby("land"): kk, nn = int(sub.est.sum()), int(len(sub)) by_land[land] = dict(established=kk, n=nn, fraction=kk / nn, jeffreys95=jeffreys(kk, nn)) cp = d[d.CPrecord == True] g3 = cp.groupby(["Binomial", "CountryName"]).agg(est=("est", "any"), known=("known", "any"), land=("LandType", mode)) k3 = g3[g3.known] by_land_cp = {land: dict(established=int(s.est.sum()), n=int(len(s)), fraction=float(s.est.mean())) for land, s in k3.groupby("land")} estrec = d[d.est] sp = d.groupby("Binomial").agg(est=("est", "any"), known=("known", "any"), spread=("IntroMethod", lambda s: (s == "Colonised.from.introd.popn").any())) ncountry = estrec.groupby("Binomial").CountryName.nunique() spe = sp[sp.est].join(ncountry.rename("ncountry")) gc = d.groupby(["Binomial", "CountryName"]).agg(est=("est", "any"), col=("IntroMethod", lambda s: (s == "Colonised.from.introd.popn").any())) nar = estrec.groupby(["Binomial", "CountryName"]).area.nunique().rename("nareas") e = gc[gc.est].join(nar) out = dict( source="GAVIA main data table (figshare 10.6084/m9.figshare.4234850)", n_records=int(len(d)), n_species=int(d.Binomial.nunique()), n_countries=int(d.CountryName.nunique()), records_with_unknown_status=int((d.StatusCat == "Unknown").sum()), records_with_unknown_status_fraction=float((d.StatusCat == "Unknown").mean()), status_counts=d.StatusCat.value_counts().to_dict(), species_by_country=dict(n=int(len(g)), known_outcome=int(len(k)), established=int(k.est.sum()), fraction_established_of_known=float(k.est.mean()), fraction_established_if_unknown_counted_as_failure=float(k.est.sum() / len(g))), by_land_type=by_land, by_land_type_first_introduction_records_only=by_land_cp, species_level=dict(n_species=int(len(sp)), with_known_outcome=int(sp.known.sum()), established_somewhere=int(sp.est.sum()), fraction_of_all_introduced_species=float(sp.est.sum() / len(sp)), fraction_of_species_with_known_outcome=float(sp.est.sum() / sp.known.sum())), spread_proxies=dict( established_species=int(len(spe)), with_population_recorded_as_colonised_from_introduced_population=int(spe.spread.sum()), fraction_colonised_flag=float(spe.spread.mean()), jeffreys95_colonised=jeffreys(int(spe.spread.sum()), len(spe)), fraction_established_in_2plus_countries=float((spe.ncountry >= 2).mean()), fraction_established_in_5plus_countries=float((spe.ncountry >= 5).mean()), per_country_established=int(len(e)), per_country_fraction_in_2plus_regions=float((e.nareas >= 2).mean()), per_country_fraction_in_3plus_regions=float((e.nareas >= 3).mean()), note="lower bounds on natural spread: the 'colonised' flag and region counts depend on recording effort"), ) return out def cassey(): d = pd.read_csv(fetch("cassey"), encoding="latin1") r = pd.to_numeric(d["R.Calc"], errors="coerce") nsz = pd.to_numeric(d["N.size"], errors="coerce") m = r.notna() & nsz.notna() & (nsz > 3) & (r.abs() < 1) z = np.arctanh(r[m]) w = nsz[m] - 3 zbar = float((w * z).sum() / w.sum()) # DerSimonian-Laird random effects Q = float((w * (z - zbar) ** 2).sum()) k = int(m.sum()) tau2 = max(0.0, (Q - (k - 1)) / (w.sum() - (w ** 2).sum() / w.sum())) wr = 1 / (1 / w + tau2) zre = float((wr * z).sum() / wr.sum()) se = float(np.sqrt(1 / wr.sum())) return dict(source="Cassey et al. 2018 S1 Data", n_rows=int(len(d)), n_effects_used=k, fraction_positive=float((r[m] > 0).mean()), median_r=float(r[m].median()), fixed_effect_r=float(np.tanh(zbar)), random_effects_r=float(np.tanh(zre)), random_effects_r_ci95=[float(np.tanh(zre - 1.96 * se)), float(np.tanh(zre + 1.96 * se))], published_statement="establishment success most critically depended on propagule pressures in the range of 10-100 individuals") def p3k_check(): """Record what the p3k14c database holds for the regions of this study (it turns out to be almost nothing).""" try: d = pd.read_csv(fetch("p3k14c"), low_memory=False) except Exception as e: # network failure must not stop the pipeline return dict(error=str(e)) pac = d[((d.Long < -100) | (d.Long > 150)) & (d.Lat < 30) & (d.Lat > -55)] rn = d[d.Province.astype(str).str.contains("Rapa Nui", case=False)] na = d[d.Country.astype(str).str.contains("Greenland|Iceland|Faroe", case=False)] return dict(n_dates=int(len(d)), continents=d.Continent.value_counts().to_dict(), remote_oceania_dates=int(len(rn)), remote_oceania_provinces=sorted(set(rn.Province.astype(str))), rapa_nui_youngest_oldest_14C_age_BP=[int(rn.Age.min()), int(rn.Age.max())] if len(rn) else None, north_atlantic_dates=int(len(na)), note="p3k14c has no dates from New Zealand, Hawaii (state field), the Cook, Society, Marquesas or other " "Polynesian islands apart from Rapa Nui, and none from Greenland or the Faroes; it cannot supply the chronology") def run(verbose=True): prop, tab, coefs = propagule() res = dict(propagule=prop, gavia=gavia(), cassey=cassey(), p3k14c=p3k_check()) lit = pd.read_csv(os.path.join(DATA, "curated", "invasion_literature_values.csv")) res["literature_values"] = lit.to_dict("records") # stage summary used by the synthesis res["stage_summary"] = dict( introduced_to_established=dict(birds_events=prop["overall_fraction"], birds_small_propagule_2_10=prop["bins"][0]["fraction"], birds_100_1000=float((prop["bins"][3]["established"] + prop["bins"][4]["established"]) / (prop["bins"][3]["n"] + prop["bins"][4]["n"])), gavia_species_lower=res["gavia"]["species_level"]["fraction_of_all_introduced_species"], gavia_known_upper=res["gavia"]["species_by_country"]["fraction_established_of_known"], literature_range=[0.10, 0.50]), established_to_spreading=dict(gavia_colonised_flag_lower=res["gavia"]["spread_proxies"]["fraction_colonised_flag"], gavia_2plus_regions=res["gavia"]["spread_proxies"]["per_country_fraction_in_2plus_regions"], literature_range=[0.10, 0.50]), ) dump_json(res, "invasion_results.json") if verbose: print("Propagule data: %d events, %.1f%% established; slope %.2f +/- %.2f per dex" % ( prop["n_events"], 100 * prop["overall_fraction"], prop["logistic_all"]["slope_per_dex"], prop["logistic_all"]["se_slope"])) print(tab[["bin", "n", "established", "fraction"]].to_string(index=False)) g = res["gavia"] print("GAVIA: species x country established %d of %d known (%.2f); by land type %s" % ( g["species_by_country"]["established"], g["species_by_country"]["known_outcome"], g["species_by_country"]["fraction_established_of_known"], {k: round(v["fraction"], 3) for k, v in g["by_land_type"].items()})) print(" spread proxies:", {k: (round(v, 3) if isinstance(v, float) else v) for k, v in g["spread_proxies"].items() if k != "note"}) print("Cassey meta-analysis: random-effects r = %.2f, %d effects, %.0f%% positive" % ( res["cassey"]["random_effects_r"], res["cassey"]["n_effects_used"], 100 * res["cassey"]["fraction_positive"])) print("p3k14c:", res["p3k14c"].get("remote_oceania_dates"), "Remote Oceania dates,", res["p3k14c"].get("north_atlantic_dates"), "North Atlantic") return res, tab, coefs if __name__ == "__main__": run()