Using Statistical Methods for Environmental Science and Management

1. Using Statistical Methods for Environmental Science and Management Graham McBride, NIWA, Hamilton g.mcbride@niwa.co.nz Statistics Teachers’ Day, 25 November 2008 What do statisticians really do?

2. THE ROLE OF STATISTICAL METHODS: MY VIEW • Separate randomness from pattern • Make inferences about the world, based on data from samples • Help to design sampling programmes (use resources efficiently) • Help to establish cause and effect • Can’t “prove anything with statistics”

3. “Three kinds of lies” Insult, or compliment? There are three kinds of lies • lies, damned lies, and statistics Who said that? • Mark Twain (1835 – 1910) “Figures often beguile me, particularly when I have the arranging of them myself” • Benjamin Disraeli (1804 – 1881) Sought to discredit true British soldier casualty figures in the Crimean War (1853 – 1856) Who came first? (Twain cites Disraeli!)

4. What you should do • Establish the context of your work (what do people want to know, and why do they want to know that?) • Consult with others, e.g., to discuss whether a proposed sampling programme can actually be done • Discuss the appropriate burden-of-proof (e.g., drinking water standards minimise the consumer’s risk, not the producer’s risk)

5. What you should not do • Confuse association and causation (pp. 267-8 of Barton, Sigma Mathematics) • Ignore other lines-of-evidence (Bradford-Hill criteria), such as • Can the cause reach the location of the effect? • Is the finding plausible? • Can you explain inconsistencies with other evidence? • Be ignorant of how statistical procedures work • The computer said so

6. What you should not do • Believe that there is only one “statistically correct” way of analysing data • There are lots of good ways; many more bad and wrong ways too • Not consider bias and imprecision in your data

7. Bias and Imprecision

8. What you might have to do • Use non-standard methods, e.g., • non-parametric (rank) methods for highly skewed data (very common in aquatic studies) • e.g., linear trend or monotonic trend? • Read rather widely • Statistics is not a cut-and-dried subject; there are still some fundamental debates about statistical inference, especially the Bayesians versus the frequentists—both approaches have their place

9. What you also might have to do • Answer this question: “What is P” • Result of a hypothesis test • Used (over-used!) routinely, so you’ll need to know • P = Prob(data at least as extreme if the tested hypothesis is true) • Not the probability of the truth of the hypothesis • Relate results to confidence intervals

10. EXAMPLEIncreasing pressure on freshwaters Is there evidence of associated deterioration (or improvements) in rivers?

11. 600000 4 Total Nitrogen 3.5 Total Phosphorus 500000 Cows 3 400000 2.5 Cow numbers (millions)2 Fertilizer consumption (tonnes)1 300000 2 1.5 200000 1 100000 0.5 0 0 1988 1990 1992 1993 1995 1996 1998 1999 2001 2002 1989 1991 1994 1997 2000 2003 Data source: 1Fertilizer consumption – UN Food & Agriculture Organisation 2Cows –Livestock Improvement NZ Dairy Statistics

12. A National River Water Quality Network for New Zealand (1989) GOAL To provide scientifically defensible information on the important physical, chemical, and biological characteristics of a selection of the nation’s rivers as a basis for advising the Minister of Science and other Ministers of the Crown of the trends and status of these waters OBJECTIVES • Detect significant trends in water quality • Develop better understanding of water resources, and hence to better assist their management

13. NRWQNstructure • 77 sites on 35 rivers • All sites have reliable flow data • Sites are sampled by regional Field Teams • 14 WQ parameters (monthly) • Data available (search for WQIS www.niwa.co.nz

15. WQ state & land use Correlations with % Pasture Temperature 0.50*** Conductivity 0.55*** pH -0.19 Dissolved oxygen -0.17 Visual clarity -0.60*** NOx-N 0.71*** NH4-N 0.77*** Total nitrogen 0.84*** DRP 0.67*** Total phosphorus 0.74*** E. coli 0.79*** ***P < 0.001; Spearman rank correlation

16. WQ Trends 1989-2005 • Calculated annual medians from monthly data at each site for each parameter • Took the 77 datapoints for each year and calculated the 5th, 50th, and 95th percentile values • The 50th percentile gives us a picture of what is happening in a national “average” river in terms of annual median water quality data • The 5th and 95th percentiles tell us about changes over time in our “best” and “worst” rivers. • Trends in these values were assessed using the Spearman rank correlation coefficient (rS).

17. 5th 50th 95th NOx-N Trends 1989-2005 1200 1000 ) 3 800 -N (mg/m 600 x NO 400 200 0 1993 1994 1995 1996 1997 1998 1999 2000 1989 1990 1991 1992 2001 2005 1988 2002 2003 2004 Year Concentrations of NOx-N increased dramatically between 1989 & 2005 in our most enriched rivers

18. Trends 1989-2005 • Results indicative of: • Warming in our coolest rivers • Drops in pH • Increasing nitrogen enrichment • Decreases in BOD5 most rivers

19. Trends 1989-2003 • More formal analysis of trends carried out on monthly data (1989-2003) at all 77 sites • Seasonal Kendall test • Data were flow-adjusted using LOWESS (many WQ parameters can be strongly influenced by discharge) • Used a binomial test to indicate a “national trend” • Discriminate between “significant” (i.e. P < 0.05) and “meaningful” trends (i.e., P < 0.05 and slope > 1% of median value per annum).

20. Trends in TN Total nitrogen exhibited a strong increasing trend at the national scale during 1989-2003 (P < 0.001). Increasing trends in TN were particularly evident in the South Island, where 25 of 33 sites showed meaningful increases.

21. Trends in DRP There was a strong national trend of increasing DRP concentrations during 1989-2003 (P < 0.001). This result contrasts with the relatively weak trends observed for 1989-2005.

22. 15 10 5 Plot 1 RSKSE 0 -5 -10 -15 Temp Cond pH DO Clar NOx-N NH4-N TN DRP TP BOD5 Summary of trends 1989-2003 No significant trend Significant improving trend Significant deteriorating trend

23. 10 y = 0.0406x - 0.0027 8 2 R = 0.31 6 4 as % of median) Trend in Dissolved Reactive Phospohorus (SKSE 2 0 -2 Lower Manawatu Rv. -4 0 10 20 30 40 50 60 70 80 90 100 % Pastoral land use Links between land use and trends The magnitude of trends in DRP increase with % pastoral land use

24. Land use and trends Parameter SKSE RSKSE Temperature 0.19 0.20 Conductivity 0.47 0.40 pH -0.28 -0.28 Dissolved oxygen -0.27 -0.27 Visual clarity -0.26 -0.11 Oxidised nitrogen 0.30 0.23 Ammoniacal nitrogen 0.29 0.68 Total nitrogen 0.35 -0.01 Dissolved reactive phosphorus 0.59 0.48 Total phosphorus 0.31 0.18 Spearman rank correlation coefficients (bold P < 0.01)

25. Conclusions • Strong associations between nutrient concentrations and %pastoral land cover at the national scale (State) • Rivers draining large areas of pastoral land have deteriorated significantly over the last 17 years with respect to nitrogen concentrations (Trends) • The magnitude of trends in some parameters is associated with extent of pastoral land use • Decreasing trends in NH4-N and BOD5 indicative of improvements in point source management • Increasing trends in nutrients indicative of increasing pressure from agriculture

26. EXAMPLE:Water quality-human health risk assessment, quantitative approachChristchurch City Wastewater Outfall

28. Quantitative Microbial Health Risk Assessment (QMHRA) • Identify hazards (pathogens) • Quantify exposure (swimming, shellfish consumption) • Assess dose-response • Characterise risk

29. Hazard vs. Risk • Hazards can cause harm, after exposure • Risk cannot occur if no exposure • Can have hazard without risk • But not vice versa!

30. Christchurch hazards—viruses only From an extensive list (next slide): • Swimming • adenovirus (respiratory) • rotavirus • enterovirus (Echovirus 12) • Shellfish consumption (raw) • enteroviruses • rotavirus • hepatitis A

32. Dose-response curves

33. Accounting for variability and uncertainty • Exposure is variable • e.g., individuals’ swim duration • Dose-response is uncertain • only some pathogen strains in clinical trials • trials limited to healthy adults • Describe using statistical distributions in a Monte Carlo analysis

34. Scenariosis! • 1,000 people; 1,000 occasions • 8 beaches • 2 influent virus conditions (normal & outbreak) • 2 seasons summer/winter • 3 viruses for 2 activities • 2 outfall lengths • 2 virus inactivation regimes • 2 UV options (with & without) •  1536 x 106 calculations

35. Calculation sequence

36. - rd = - Pr ( d ) 1 e inf Dose-response models • Constant susceptibility—simple exponential (d = average dose, Prinf = infection prob) • Variable susceptibility—“beta-Poisson” • Calculations performed using “@RISK” (an Excel plug-in)

37. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 1, Individual 1 Volume ingested Dose Probability of infection Binomial distribution Infected?

38. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 1, Individual 2 Volume ingested Dose Probability of infection Binomial distribution Infected?

39. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 1, Individual 3 Volume ingested Dose Probability of infection Binomial distribution Infected?

40. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Sum the cases Occasion 1, Individual 1000 Volume ingested Dose Probability of infection Binomial distribution Infected?

41. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 2, Individual 1 Volume ingested Dose Probability of infection Binomial distribution Infected?

42. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 2, Individual 2 Volume ingested Dose Probability of infection Binomial distribution Infected?

43. Duration Frequency Prob(inf) Ingestion rate Microorg. concn Frequency Frequency Occasion 2, Individual 3 Volume ingested Dose Probability of infection Binomial distribution Infected?

44. Characterising the results • Risk percentiles—percent of time the risk is below a stated value • IIR—Individual Infection Risk (total number of calculated infections divided by total number of exposures)

45. Results South New Brighton Integers are cases per 1000 exposures

46. IIR: Normal influent, South Brightonadenovirus, swim Numbers are percentages. MfE/MoH (2003) guidelines: <0.3% = “Very good”.

47. IIR: Normal influent, South Brighton rotavirus, shellfish Numbers are percentages.

48. IIR: Outbreak influent, South Brighton adenovirus, swim Numbers are percentages.MfE/MoH (2003) guidelines: 1.9 - 3.9% = “Fair” - “Poor”.

49. IIR: Outbreak influent, South Brighton rotavirus, shellfish Numbers are percentages.

50. IIR: Outbreak influent, South Brighton hepatitis A, shellfish Numbers are percentages.