Mixed-model analysis of a percentage outcome using logit transformation in R

Warning: this post is terribly technical. Readers not inclined towards statistical analysis and R coding may wish to skip this one!

I was recently contacted by a biologist in the UK regarding one of the analyses that I did in my PhD. This researcher was trying to model percentage data with mixed effects models and was struggling to find examples of where this had been done before to base his analysis on, except for the analysis that I performed to identify drivers of disturbance at a patch level in the Great Western Woodlands. I remember that at the time, I really struggled to figure out the best way to perform the analysis, it was so complicated and there were issues that made it different from all the other examples I was trying to go off. I remember assuming that it was probably really straightforward and that everyone else does these sorts of analyses all the time and that I was a bit slow for not being able to figure it out straight away… so it is somewhat gratifying to hear that it’s not a common problem that biologists solve, and to have someone else reaching out to me for my experience on it!

I’ve decided to write up my response as a blog post, in case anyone else is interested… and to remind myself of what happened there! Comments and suggestions for improvement are invited!

In general terms, the analysis is aimed at understanding what factors explain how much anthropogenic disturbance (e.g. land clearing) will occur within the landscapes of the Great Western Woodlands, where there is an abundance of mining activity and pastoralism within an otherwise relatively undisturbed region.

More specifically, the analysis relates to predicting the percentage of area disturbed (the ‘development footprint’) at a patch level (this was the smaller of two scales of analysis used in the investigation, and the only one to be discussed here), by a set of potential predictor variables. The predictor variables included both categorical and continuous variables. Categorical variables were: pastoral tenure, greenstone lithology, conservation estate, ironstone formation, schedule 1 area clearing restrictions, environmentally sensitive area designation, vegetation formation, type of tenements (exploration/prospecting tenement, mining and related activities tenement, none), primary target commodity (gold, nickel, iron-ore, other), and sample area (n=24). The numerical variables were: number of mining projects, number of dead tenements, sum of duration of all live and dead tenements, distance to wheatbelt, and distance to nearest town.

The next image (below) shows an example of these patches and their attribuets.

pic for blog 2.jpg

Example of polygons constructed to assess predictors of development footprint at the patch level

I used ‘nlme’ package in R for the analysis. Early in the data exploration (which I won’t detail here) I realised that the data for the proportional area disturbed was very skewed, and that log or other transformations weren’t enough to achieve a normal-ish distribution. A logistic transformation was the clear solution, given the percentage nature of the data (i.e. that it was bounded by 0 and 1); see the following image with scatter plots and histograms.

pics-for-blog-post

However, I couldn’t just use a GLMM with a gaussian distribution and a logit link, as this would have produced the problem that predicated values will be outside the 0-1 range. I had to first do a logistic transformation, model linearly, and then back-transform the modelled values.

The following line of code shows how I logit-transformed the proportion of area disturbed (‘dist_perc’), after adding 0.001 to each value (I had to do this as logit only works on values between 0 and 1, exclusive). I also show how I ‘untranformed’ the data and plotted it to produce a nice, normal, histogram. Note that dp was the dataframe I worked with, stands for ‘disturbance by patch’.

> dp$ldistp = qlogis(dist_perc+0.001)                   # qlogis does a logit transformation on (a set of) values.

> plot(dp$ldistp)                                         # scatter plot looks much better

> range(dp$ldistp)

> hist(dp$ldistp)                                         # Looks good.

> dp$s.dist = scale(dp$ldistp)                   # s.dist is the logit-transformed, scaled disturbed percentage

> back_transformed = (plogis(dp$ldistp))-0.001           # plogis ‘undoes’ the logit transformation.

 

After completion of the data exploration and dealing with other variables that needed to be transformed, confirming no outliers, dealing with collinearity etc and scaling issues, I modelled the transformed variable using maximum likelihood and performed a manual backward-selection and a few other checks to reach the final model, which I then modelled using REML.

I have included the whole script here. It includes the date exploration stages, as well as a plethora of comments, and also a few tries of alternative analysis techniques using lmer, lmerTest and lme4 (and modelling with a logit-link before I figured out that this wasn’t a useful approach). I take no responsibility for errors! Check line 621 for the final model, followed by back-transformations of the outcomes.

Advertisements
Posted in blog | Tagged , , | 1 Comment

R Script for spatial analysis disturbance by patch – for reference in blog post on mixed-model analysis with logistic transformation

Note: this is a post to hold the R script that is referred to in this post. It’s even more awfully technical than the post that refers to it; readers not interested or familiar with R code or this type of analysis may want to skip it!

You’d probably want to copy this all in to an R script (and perhaps view in r studio) to make sense of it.

# Disturbance by patch to determine landscape predictors of disturbance —-

# Housekeeping —-
setwd(“C:/Users/ke/Dropbox/PhD/Chapter 3 – spatial analysis/R”)
library(lme4); library(ggplot2); library(plyr); library(lattice); library(car); require(lmerTest); library(nlme)
source(“HighstatLibV6.R”)
dp

# Start of data exploration —-
names(dp)
str(dp)
summary(dp)
attach(dp)
head(dp)
sum(is.na(dp)) # no missing values

# 4.2 looking for outliers
vars1 “greenstone” , “esa” , “cons_estate” , “iron_formation”)
Mydotplot(dp[, vars1])
vars2 “mining_activity”,”commodity”,”wheatbelt_dist”,”town_dist”,”patch_area”,”disturbed.area”, “dist_perc”)
Mydotplot(dp[, vars2])
vars3 Mydotplot(dp[, vars3])
# no outliers once the tenement years issue was fixed (in raw data), but some vars need transformation

# Transforming variables —-

# mini-tutorial on how scaling works:
x s.x s.x # the resultant scaled variable stored the scaling attributes: the centre (mean) and scale (multiplier)
attr(s.x,”scaled:center”) # 5.5; extracting the centre, the mean of the original vector, from which each value was subtracted to do the scaling
attr(s.x,”scaled:scale”) # 3.02765; extracting the scale (the multiplier)
x.recreated = s.x*attr(s.x,’scaled:scale’) + attr(s.x,’scaled:center’) # to unscale just multiply the scaled vector by the ‘scale’ and add the ‘centre’

# transforming disturbance percentage:
plot(dp$polyID~dp$dist_perc)
hist(dp$dist_perc) # too skewed
dist_perc1 = dist_perc + 1
dp$logdisp = log(dist_perc1)
range(dp$logdisp)
hist(dp$logdisp)
plot(dp$polyID~dp$logdisp) # still too squashed. try log10 transformation?
dp$logdisp = log10(dist_perc1)
attach(dp)
range(dp$logdisp)
hist(dp$logdisp)
plot(dp$polyID~dp$logdisp) # still too skewed. Let’s try a logit transformation:
qlogis(0.9) # Getting familiar w logit transformations… this returns 2.197. qlogis does a logit transformation on (a set of) values.
plogis(2.197225) # returns 0.9. plogis ‘undoes’ the logit transformation.
dp$ldistp = (qlogis(dp$dist_perc+0.001)) # ldistp is logit transformation of the percent disturbed plus 0.001
plot(dp$ldistp) # Looks good.
range(dp$ldistp)
hist(dp$ldistp)
dp$s.dist = scale(dp$ldistp) # s.dist is the logit-transformed, scaled disturbed percentage
hist(dp$s.dist)
# OK I need to either do a logistic regression here, or simply logit-transform the response before modelling linearly
# doing a logistic transformation and then modelling linearly – advised by Michael Renton

# transforming and scaling ten_yrs:
range(ten_yrs) # the 177 is a valid value.
plot(ten_yrs) # needs some transformation
dp$lntyrs = log((ten_yrs+1))
plot(dp$lntyrs) # yep the natural log looks good.
range(dp$lntyrs)
dp$s.lntyrs = scale(dp$lntyrs, center = TRUE, scale = TRUE)
plot(dp$lntyrs)
plot(dp$s.lntyrs)

# transforming patch area:
range(dp$patch_area)
dp$patchkm2 = dp$patch_area/1000000 # make kms instead of m2
range(dp$patchkm2)
hist(dp$patchkm2)
dp$ln.pkm = log(dp$patchkm2) # natural log of patch area
hist(dp$ln.pkm) # centred and normalish

# transforming minedex_count:
plot(minedex_count)
table(dp$minedex_count)
plot(log(minedex_count+1)) # somewhat better
logmx = log(minedex_count+1)
range(logmx)
dp$s.mx = scale(logmx)
plot(dp$s.mx)

# transforming distance to wheatbelt:
range(wheatbelt_dist)
dp$wbkm = wheatbelt_dist/1000
plot(dp$wbkm)
range(dp$wbkm)
dp$s.wb = scale(dp$wbkm)
plot(dp$s.wb)

# transforming distance to town:
dp$tnkm = town_dist/1000
plot(dp$tnkm)
range(dp$tnkm)
dp$s.tn = scale(dp$tnkm)
plot(dp$s.tn)

# not much can be done with the commodity variable:
summary(dp$commodity) # not sure if having the “no commodity stated” value is useful…
# dp$commodity[dp$commodity==”no commodity stated”] <- “none” # creates NAs

# adding second-order polynomials of the numerical variables
dp$tenyrs2 = dp$s.lntyrs^2
dp$mx2 = dp$s.mx^2
dp$wb2 = dp$s.wb^2
dp$tn2 = dp$s.tn^2
range(dp$s.lntyrs)
range(dp$tenyrs2)

# combining greenstone and ironstone lithology data into one categorical variable
dp$lithology = ifelse(greenstone == “greenstone”,”greenstone”,”other.lithology”)
dp$iron_formation
dp$lithology = ifelse(iron_formation == “iron-form”,”iron.formation”,dp$lithology)
dp$lithology
table(dp$lithology) # clever duck

# Looking for collinearity —-
pairs(dp[,vars1], lower.panel = panel.cor)
#i.e. SA number, SA_mines and category all related – use SA_mines
pairs(dp[,vars2], lower.panel = panel.cor)
miningvars = c(“SA_mines”,”greenstone”,”iron_formation”,”minedex_count”,”ten_yrs”,”mining_activity”,”commodity”)
pairs(dp[,miningvars], lower.panel = panel.cor)
# also strong collinearity between minedex count and dt-count, mined and
# note mining activity replaces expl_ten and mined.
# low 0.4 collinearity between minedex count and ten_yrs, and 0.5 between SA_mines and mining activity

summary(dp$dist_perc ~ dp$veg_form)
summary(dp$cons_estate)

dpvars = c(“pastoral”,”greenstone”,”esa”,”cons_estate”,”iron_formation”,”sched_1″,”veg_form”,
“mining_activity”,”commodity”,”s.mx”,”s.lntyrs”,”s.wb”,”s.tn”,”s.dist”)
Mydotplot(dp[, dpvars])
pairs(dp[,dpvars], lower.panel = panel.cor)
# ok, have dealt with collinearity. There’s nothing above 0.5, and even the 0.5 is not important

# looking at possible X Y relationships —-

# first the numeric variables: plot Y against the Xes and calculate pearson coefficients:
dpnumvars = c(“s.mx”,”s.lntyrs”,”s.wb”,”s.tn”,”s.dist”)
pairs(dp[,dpnumvars], lower.panel = panel.cor)
# there appears to be a positive relationship between the area disturbed and both the
# number of mining projects and the total duration of tenements in an area.
# No relation with distance from wheatbelt, and slight negative relationship with distance from town,
# which makes sense – areas near towns are more likely to get more disturbed and to stop tracks from growing over etc.

plot(dp$dist_perc~dp$lntyrs) #remember that dp$lntyrs = log((ten_yrs+1))
plot(dp$ldistp~dp$lntyrs)

#then categorical variables:
boxplot(dp$s.dist~dp$pastoral) # no real effect of pastoral activity
boxplot(dp$s.dist~dp$greenstone) # more disturbance on greenstone
boxplot(dp$s.dist~dp$esa) # more disturbance on ESAs
boxplot(dp$s.dist~dp$cons_estate) # greatest disturbance on former leasehold then class As!
# UCL has same mean distrubance as gazetted conservation estate; though higher variance.
boxplot(dp$s.dist~dp$iron_formation) # more disturbance on BIFs
boxplot(dp$s.dist~dp$sched_1) # more disturbance on schedule 1 areas
boxplot(dp$s.dist~dp$veg_form) # less disturbance in bare areas and succulent steppe
boxplot(dp$dist_perc~dp$veg_form)
boxplot(dp$s.dist~dp$mining_activity) # more disturbance where there’s been mining tenement.
# strangely similar between exploration and no mining tenement.
boxplot(dp$s.dist~dp$commodity) # most disturbance around gold mining projects, maybe not significant

# Looking for interactions —-

coplot(s.dist ~ s.mx | factor(commodity), data = dp)
# R displays from lower left, across, to upper right.

coplot(s.dist ~ s.mx | factor(commodity), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# there could be an interaction between the commodity and the number of mining projects

coplot(s.dist ~ s.lntyrs | factor(pastoral), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# There could be an interaction between tenement years and pastoral, same as found in regional-scale analysis.

coplot(s.dist ~ factor(greenstone) | factor(pastoral), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# maybe nothing between greenstone and pastoral

coplot(s.dist ~ factor(veg_form) | factor(greenstone), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# maybe nothing between greenstone and pastoral
coplot(s.dist ~ s.tn | factor(veg_form), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# well there doesn’t seem to be much mallee near towns; looks like there could be an interaction

coplot(s.dist ~ s.tn | factor(esa), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# definitely looks like an interaction between esas and dist to town

coplot(s.dist ~ s.wb | factor(esa), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })
# same again – interaction between esas and distance to wheatbelt, although esas tend to be much closer to wheatbelt than non-esas.

coplot(s.dist ~ factor(veg_form) | factor(esa), data = dp,
panel = function(x, y, …) {
tmp abline(tmp)
points(x, y) })

coplot(s.dist ~ factor(mining_activity) | factor(veg_form), data = dp )
# most of the mining and exploration appears to be in woodland

# Fitting and selecting models —-

# 4.7.1 using gaussian family and logit link with lme4
# glmer example: model1 dp1 = glmer(dist_perc ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+
(1|SA_number),
data=dp, family = gaussian (link = ‘logit’))
summary(dp1)

E1 F1 plot(x = F1,
y = E1,
xlab = “Fitted values”,
ylab = “Residuals”,
main = “Homogeneity?”)
abline(h = 0, v = 0, lty = 2)
plot(dp1)
plot(F1)

# 4.7.2 Fitting a linear mixed model with a logit-transformed response in lmerTest
dp2 = lmer(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+tn2+
(1|SA_number), REML=TRUE,
data=dp)
summary(dp2)
plot(dp2)
summary(dp$commodity)
# dp$commodity[is.na(dp$commodity)] <- “none”

# compare with model that doesn’t have random effects to see if they are important
dp2.gls = gls(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+tn2,
data=dp)
summary(dp2.gls)
plot(dp2.gls)

anova(dp2.gls,dp2)
anova(dp2)

plot(dp2)
E2 F2 plot(x = F2,
y = E2,
xlab = “Fitted values”,
ylab = “Residuals”,
main = “Homogeneity?”)
abline(h = 0, v = 0, lty = 2)
table(dp$SA_number)
length(E2)

dp3 = lmer(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+tn2+
(1|SA_number),
data=dp, REML = FALSE)
summary(dp3)
anova(dp3)

dp4 = lmer(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+
(1|SA_number),
data=dp, REML = FALSE)
summary(dp4)
anova(dp4)

dp5 = lmer(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+
(1|SA_number),
data=dp, REML = FALSE)
summary(dp5)
anova(dp5)

dp6 = lmer(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+
(1|SA_number),
data=dp, REML = FALSE)
summary(dp6)
anova(dp6)

# 4.7.3 Fitting a linear mixed model with a logit-transformed response in nlme

library(nlme)
# first check the significance of the random effect structure:
dp1 = lme(s.dist ~ factor(pastoral)+factor(iron_formation)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(sched_1)+factor(veg_form) +factor(commodity)+factor(mining_activity)+
s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+tn2+ln.pkm,
random = ~1|SA_number,method = “REML”,
data=dp)
summary(dp1)
plot(dp1)

# compare with model that doesn’t have random effects to see if they are important
dp1.gls = gls(s.dist ~ factor(pastoral) +factor(iron_formation) +factor(greenstone) +factor(esa) +
factor(cons_estate) +factor(sched_1) +factor(veg_form) +factor(commodity) +
factor(mining_activity) +s.lntyrs +s.mx +s.wb +s.tn +tenyrs2 +mx2 +wb2 +tn2 +ln.pkm,
data=dp)
dp$s.dist
summary(dp1.gls)
plot(dp1.gls)

anova(dp1.gls,dpn1)
# The p-value of the Maximum Likelihood ratio test comparison between the mixed and fixed methods
# (both calculated using REML, which allows us to apply the likelihood ratio test) is 0.0154,
# indicating that the random factor is significant.

# check the model:
E2 F2 op MyYlab <- “Residuals”
plot(x = F2, y = E2, xlab = “Fitted values”, ylab = MyYlab, main=”Residuals versus fitted”)
boxplot(E2 ~ pastoral, data = dp,
main = “Pastoral tenure”, ylab = MyYlab)
boxplot(E2 ~ greenstone, data = dp,
main = “Greenstone lithology”, ylab = MyYlab)
plot(x = dp$s.lntyrs, y = E2, ylab = MyYlab,
main = “Tenement duration”, xlab = “scaled log of total tenement duration”)
par(op)
#haven’t checked against all variables but there’s too many, moving on for now and will check when model is smaller.

# Next, I’m going to fit the ‘beyond optimal’ model with ML, then reduce it manually
# then validate it,
dp2 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+tn2+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp2)

# remove tn2
dp3 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+mx2+wb2+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
anova(dp2,dp3)

# remove mx2
dp4 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+wb2+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
anova(dp4,dp3) #yes, wasn’t significant
summary(dp4)

# remove wb2
dp5 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+tenyrs2+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp5)

#remove tenyrs2
dp6 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.mx+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp6)

# remove s.mx
dp7 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(sched_1)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp7)

# remove schedule 1 areas
dp8 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(iron_formation)+factor(veg_form)+factor(commodity)+
factor(mining_activity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp8)

# remove iron_form
dp9 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+factor(mining_activity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp9)

# remove either distance to wheatbelt or mining activity
dp9.no.wb = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+factor(mining_activity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
dp9.no.ma = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
anova(dp9,dp9.no.wb) # p = 0.1054
anova(dp9,dp9.no.ma) # p = 0.2048

# so remove mining activity

dp10 = lme(s.dist ~ factor(pastoral)+factor(greenstone)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp10)

# not significant are: greenstone, cons estate, wheatbelt. try remove each, one at a time.
dp10.no.green dp10.no.cons dp10.no.wb

anova(dp10,dp10.no.green) # p = 0.0539
anova(dp10,dp10.no.cons) # p < 0.0001
anova(dp10,dp10.no.wb) # p = 0.0395

# got to remove greenstone
dp11 = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.wb+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp11)

# pastoral is significant
# esa is significant
# cons estate need to be checked but is probably highly significant
# veg form looks significant
# commodity looks significant
# tenement years is highly significant
# s.wb – what is that still doing here? the p values from the model are so different from anova comparing models!

dp11.no.wb = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
anova(dp11,dp11.no.wb)

# yep get rid of wheatbelt
dp12 = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp12)

# I think we’re there so just double checking the significance of every variable:

dp12.no.pastoral = lme(s.dist ~ factor(esa)+factor(cons_estate)+factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.esa = lme(s.dist ~ factor(pastoral)+factor(cons_estate)+factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.cons = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.veg = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.comm = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.lntyrs = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+factor(commodity)+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.tn = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+factor(commodity)+s.lntyrs+ln.pkm,
random = ~1|SA_number,method = “ML”, data=dp)
dp12.no.pkm = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+factor(commodity)+s.lntyrs+s.tn,
random = ~1|SA_number,method = “ML”, data=dp)

anova(dp12,dp12.no.pastoral) # p = 0.0046
anova(dp12,dp12.no.esa) # p = 0.0001
anova(dp12,dp12.no.cons) # p = 0.0001
anova(dp12,dp12.no.veg) # p = 0.0206
anova(dp12,dp12.no.comm) # p = 0.0001
anova(dp12,dp12.no.lntyrs) # p = 0.0001
anova(dp12,dp12.no.tn) # p = 0.000001
anova(dp12,dp12.no.pkm) # p = 0.0348
# all significant, and lower AIC values!

# what if i try to reintroduce greenstone?
dp12.plus.green = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(greenstone)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
anova(dp12,dp12.plus.green)
# greenstone is marginally significant with p = 0.0562 and a slightly lower AIC value (but within 2 AIC). Leave out.

#what about interaction between mining and pastoralism?
dp13 = lme(s.dist ~ factor(pastoral)+factor(esa)+factor(cons_estate)+
factor(veg_form)+factor(commodity)+s.lntyrs+s.tn+ln.pkm+s.lntyrs:factor(pastoral),
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp13)
# interaction not significant

# I’ll also just try to use non-scaled versions of the variables for ease of interpretation if possible:
dp14 = lme(ldistp ~ factor(pastoral)+factor(esa)+factor(cons_estate)+
factor(veg_form) +factor(commodity) +lntyrs +tnkm +ln.pkm,
random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp14)
# excellent! unscaled it is…

# what about removing commodity as it’s not helpful as is, and seeing if this
# and allows for reinstatement of greenstone too, and removal of patch area:
dp15 = lme(ldistp ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+
lntyrs+tnkm+factor(greenstone),random = ~1|SA_number,method = “ML”,
data=dp)
summary(dp15)
#yes, lets just recheck the significance of all variables:

dp15.no.pastoral = update(dp15, .~. – factor(pastoral))
dp15.no.esa = update(dp15, .~. – factor(esa))
dp15.no.cons = update(dp15, .~. – factor(cons_estate))
dp15.no.veg = update(dp15, .~. – factor(veg_form))
dp15.no.green = update(dp15, .~. – factor(greenstone))
dp15.no.lntyrs = update(dp15, .~. – lntyrs)
dp15.no.tn = update(dp15, .~. – tnkm)
dp15.w.pkm = lme(ldistp ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+lntyrs+
tnkm+ln.pkm+factor(greenstone)+ln.pkm,random = ~1|SA_number,method = “ML”,data=dp)

anova(dp15,dp15.no.pastoral) # p = 0.0031
anova(dp15,dp15.no.esa) # p = 0.0004
anova(dp15,dp15.no.cons) # p < 0.0001
anova(dp15,dp15.no.veg) # p = 0.005
anova(dp15,dp15.no.green) # p = 0.006
anova(dp15,dp15.no.lntyrs) # p = 0.0001
anova(dp15,dp15.no.tn) # p = 0.0002
anova(dp15,dp15.w.pkm) # p = 0.778 # excluding patch area is justified.
# otherwise all significant, and lower AIC values!

# Final dp model, using REML and with unscaled response variable
finaldp = lme(ldistp ~ factor(pastoral)+factor(esa)+factor(cons_estate)+factor(veg_form)+
factor(greenstone)+lntyrs+tnkm,random = ~1|SA_number,method = “REML”, data=dp)
summary(finaldp)

# Validate model
par(mfrow=c(2,2))
plot(finaldp)
E1 F1

par(mfrow = c(1, 1))
plot(x = F1,
y = E1,
xlab = “Fitted values”,
ylab = “Residuals”,
main = “Homogeneity?”)
abline(h = 0, v = 0, lty = 2)
# according to Zuur 2009, looks ok
boxplot(E1 ~ factor(dp$greenstone)) # compare variance across factor levels
boxplot(E1 ~ factor(dp$esa))
boxplot(E1 ~ factor(dp$cons_estate))
boxplot(E1 ~ factor(dp$veg_form))
boxplot(E1 ~ factor(dp$commodity))
boxplot(E1 ~ factor(dp$mining_activity))
tapply(E1, INDEX = dp$pastoral, FUN = var)
tapply(E1, INDEX = dp$cons_estate, FUN = var)
# variance varies only by a factor of ~4; Zuur says that it’s a problem if greater than 10. So not a problem.
# could consider removing the outlier

# test for independence:
plot(x = dp$patch_area,
y = E1,
xlab = “patch area”,
ylab = “Residuals”)
abline(h = 0, lty = 2)

cor.test(E1, dp$patch_area) # not significant

# i.e. no clusters above or below; looks ok.

#Normality

hist(E1, main = “Normality”, breaks=10)
#Or qq-plot
qqnorm(E1)
qqline(E1)
# i.e. not so normal
hist(dp$s.dist ) # very skewed no worries

#But plot residuals also versus covariates NOT in the model!
boxplot(E1 ~ s.wb) ; abline(h = 0, lty = 2) # confirms no effect of distance to wheatbelt
boxplot(E1 ~ s.tn) ; abline(h = 0, lty = 2) # no residual effect of distance to town

# Final model and sketch with predictions —-

finaldp = lme(ldistp ~ factor(pastoral) +factor(esa) +factor(cons_estate) +factor(veg_form) +
factor(greenstone)+lntyrs+tnkm,random = ~1|SA_number,method = “REML”, data=dp)
summary(finaldp)
summary(finaldp)$coef
anova.lme(finaldp)

# I tried many different ways to predict from the model using a newdata dataframe but couldn’t get it to work.
# here i use the values i calculated manually in excel using the model outputs:

# remember how to un-logit the response variable:
plogis(2.197225) # returns 0.9. plogis undoes the logit transformation.
ldistp = (qlogis(dist_perc+0.001))

# predicting effect of pastoral status:
pastoral.pred.logit = data.frame(pastoral = c(“pastoral”,”pastoral”,”pastoral”,”not.pastoral”,”not.pastoral”,”not.pastoral”),
logit.percent.disturbed = c(-5.48660829,
-5.29961149,
-5.67360509,
-4.94659229,
-4.75959549,
-5.13358909))
pastoral.pred.logit$percent.disturbed = plogis(pastoral.pred.logit$logit.percent.disturbed)-0.001
pastoral.pred.logit$percent.disturbed
pastoral.pred2 fit = (c(0.0031247841,0.0060574270)*100))
pastoral.pred2$se.low = (c((0.003124784-0.002423696),(0.006057427-0.004860812))*100)
pastoral.pred2$se.high = (c((0.003968722-0.003124784),(0.007496270-0.006057427))*100)
k k + geom_pointrange(size=1,aes(colour=Pastoral.status)) + theme_bw()+ylab(“Percent of area disturbed”)
ggsave(“dp figure 1 – pastoral status.png”, width=4, height=3, dpi=300)

# using intervals function to calculate predictions and 95% intervals:
conf.ints = intervals(finaldp)
# predicting effect of pastoral status:
pastoral.pred1 = data.frame(pastoral = c(“pastoral”,”pastoral”,”pastoral”,”not.pastoral”,”not.pastoral”,”not.pastoral”),
logit.percent.disturbed = c(-5.40035707,
-3.81570717,
-6.98500697,
-4.86034100,
-3.30078471,
-6.41989730))
pastoral.pred1$percent.disturbed = plogis(pastoral.pred1$logit.percent.disturbed)-0.001
attach(pastoral.pred1)
fit1=percent.disturbed[1];low1 = percent.disturbed[3];high1=percent.disturbed[2]
fit2=percent.disturbed[4];low2 = percent.disturbed[6];high2=percent.disturbed[5]
pastoral.pred2 fit = c(fit1,fit2))
pastoral.pred2$se.low = c((fit1-low1),(fit2-low2))
pastoral.pred2$se.high = c((high1-fit1),(high2-fit2))

k k + geom_pointrange(size=1,aes(colour=Pastoral.status)) + theme_bw()+ylab(“Percent of area disturbed”)
ggsave(“dp figure 1 – pastoral status(full uncertainty).png”, width=4, height=3, dpi=300)

# predicting effect of esa status:
esa.pred1 = data.frame(esa = c(“esa”,”esa”,”esa”,”not.esa”,”not.esa”,”not.esa”),
logit.percent.disturbed = c(-3.83612729,
-3.64913049,
-4.02312409,
-4.94659229,
-4.75959549,
-5.13358909))
esa.pred1$percent.disturbed = (plogis(esa.pred1$logit.percent.disturbed)-0.001)*100
attach(esa.pred1)
fit1=percent.disturbed[1];low1 = percent.disturbed[3];high1=percent.disturbed[2]
fit2=percent.disturbed[4];low2 = percent.disturbed[6];high2=percent.disturbed[5]
esa.pred2 fit = c(fit1,fit2))
esa.pred2$se.low = c((fit1-low1),(fit2-low2))
esa.pred2$se.high = c((high1-fit1),(high2-fit2))
k k + geom_pointrange(size=1,aes(colour=Clearing.regs)) + theme_bw()+ylab(“Percent of area disturbed”)
ggsave(“dp figure 2 – esa.png”, width=4, height=3, dpi=300)

# predicting effect of conservation estate:
table(dp$cons_estate)
cons.pred1 = data.frame(cons = c(rep(c(“class-A”,”ex.leasehold”,”gazetted”,”not-cons”,”unofficial”),each=3)),
logit.percent.disturbed = c(-5.70267029, -4.91130149, -6.49403909,
-4.34291529, -3.60996819, -5.07586239,
-6.59664929, -5.88199979, -7.31129879,
-4.94659229, -4.30305129, -5.59013329,
-5.45534029, -4.75924019, -6.15144039))
cons.pred1$percent.disturbed = plogis(cons.pred1$logit.percent.disturbed)-0.001
attach(cons.pred1)

fit1= percent.disturbed[1]; low1 = percent.disturbed[3]; high1=percent.disturbed[2]
fit2= percent.disturbed[4]; low2 = percent.disturbed[6]; high2=percent.disturbed[5]
fit3= percent.disturbed[7]; low3 = percent.disturbed[9]; high3=percent.disturbed[8]
fit4= percent.disturbed[10];low4 = percent.disturbed[12];high4=percent.disturbed[11]
fit5= percent.disturbed[13];low5 = percent.disturbed[15];high5=percent.disturbed[14]

cons.pred2 fit = c(fit1,fit2,fit3,fit4,fit5))
cons.pred2$se.low = c((fit1-low1),(fit2-low2),(fit3-low3),(fit4-low4),(fit5-low5))
cons.pred2$se.high = c((high1-fit1),(high2-fit2),(high3-fit3),(high4-fit4),(high5-fit5))
k k + geom_pointrange(size=1,aes(colour=Conservation.tenure)) + theme_bw()+ylab(“Percent of area disturbed”)
ggsave(“dp figure 3 – cons.png”, width=6, height=3, dpi=300) # note colouring is different to combined plot (r2 code) but content is same

# predicting effect of vegetation form:
table(dp$veg_form)
veg.pred1 = data.frame(veg = c(rep(c(“bare”,”broombush”,”mallee”,”mulga”,”shrubland”,”succulent”,”woodland”),each=3)),
logit.percent.disturbed = c(-4.56203229, -5.35340109, -3.77066349,
-4.12225829, -4.35071239, -3.89380419,
-4.33060829, -4.71315739, -3.94805919,
-3.40862529, -3.71443229, -3.10281829,
-4.10232729, -4.37808419, -3.82657039,
-4.00632829, -4.30158499, -3.71107159,
-3.88034929, -4.08146089, -3.67923769))
veg.pred1$percent.disturbed = plogis(veg.pred1$logit.percent.disturbed)-0.001
attach(veg.pred1)
fit1=percent.disturbed[1];low1 = percent.disturbed[2];high1=percent.disturbed[3]
fit2=percent.disturbed[4];low2 = percent.disturbed[5];high2=percent.disturbed[6]
fit3=percent.disturbed[7];low3 = percent.disturbed[8];high3=percent.disturbed[9]
fit4=percent.disturbed[10];low4 = percent.disturbed[11];high4=percent.disturbed[12]
fit5=percent.disturbed[13];low5 = percent.disturbed[14];high5=percent.disturbed[15]
fit6=percent.disturbed[16];low6 = percent.disturbed[17];high6=percent.disturbed[18]
fit7=percent.disturbed[19];low7 = percent.disturbed[20];high7=percent.disturbed[21]
veg.pred2 fit = c(fit1,fit2,fit3,fit4,fit5,fit6,fit7))
veg.pred2$se.low = c((fit1-low1),(fit2-low2),(fit3-low3),(fit4-low4),(fit5-low5),(fit6-low6),(fit7-low7))
veg.pred2$se.high = c((high1-fit1),(high2-fit2),(high3-fit3),(high4-fit4),(high5-fit5),(high6-fit6),(high7-fit7))
k k + geom_pointrange(size=1,aes(colour=Vegetation.form)) + theme_bw()+ylab(“Percent of area disturbed”)+xlab(“Vegetation formation”)
ggsave(“dp figure 4 – veg.png”, width=6, height=3, dpi=300)

summary(dp$dist_perc~dp$veg_form)
summary(veg.pred$percent.disturbed~ veg.pred$veg)

# predicting effect of greenstone lithology:
summary(dp$dist_perc~dp$greenstone)
green.pred1 = data.frame(green = c(rep(c(“greenstone”,”not_greenstone”),each=3)),
logit.percent.disturbed = c(-4.55671829, -5.34808709, -3.76534949,
-4.94659229, -5.09247909, -4.80070549))
green.pred1$percent.disturbed = plogis(green.pred1$logit.percent.disturbed)-0.001
attach(green.pred1)
fit1=percent.disturbed[1];low1 = percent.disturbed[2];high1=percent.disturbed[3]
fit2=percent.disturbed[4];low2 = percent.disturbed[5];high2=percent.disturbed[6]
green.pred2 fit = c(fit1,fit2))
green.pred2$se.low = c((fit1-low1),(fit2-low2))
green.pred2$se.high = c((high1-fit1),(high2-fit2))
k k + geom_pointrange(size=1,aes(colour=Lithology)) + theme_bw()+ylab(“Percent of area disturbed”)+xlab(“Lithology”)
ggsave(“dp figure 5 – greenstone.png”, width=4, height=3, dpi=300)

# predicting effect of tenement years:
range(lntyrs) # 0.000000 to 5.181784
tenyrs.pred = data.frame(lntyrs = seq(from=0, to= 5.181784, by= 0.01))
tenyrs.pred$logit.percent.disturbed = -5.107206754 + 0.440025 *tenyrs.pred$lntyrs
tenyrs.pred$logit.percent.disturbed.upper = -5.107206754 + (0.440025+0.0756683) *tenyrs.pred$lntyrs
tenyrs.pred$logit.percent.disturbed.lower = -5.107206754 + (0.440025-0.0756683) *tenyrs.pred$lntyrs
tenyrs.pred$tenyrs = exp(tenyrs.pred$lntyrs)-1
tenyrs.pred$percent.disturbed = (plogis(tenyrs.pred$logit.percent.disturbed)-0.001)*100
tenyrs.pred$percent.disturbed.upper = (plogis(tenyrs.pred$logit.percent.disturbed.upper)-0.001)*100
tenyrs.pred$percent.disturbed.lower = (plogis(tenyrs.pred$logit.percent.disturbed.lower)-0.001)*100
plot(percent.disturbed~lntyrs, data = tenyrs.pred, xlab = “log of tenement years”, ylab=”Percentage of area disturbed”)

head(tenyrs.pred) # variables of use are: tenyrs,percent.disturbed,percent.disturbed.upper,percent.disturbed.lower

#start of plot
p = ggplot(tenyrs.pred, aes(tenyrs,percent.disturbed))
p + geom_line(colour=I(“red”),size=1.1)+
geom_ribbon(aes(ymin=percent.disturbed.lower, ymax=percent.disturbed.upper),
fill=I(“red”),alpha=0.2)+
theme_bw() +
labs(x=”Total tenement duration (years)”,y=”Area disturbed (%)”)
ggsave(“dp figure 6 – tenement years.png”, width=4, height=3, dpi=300)
# end of plot

# predicting effect of distance to town:
range(tnkm) # 0.74 100.00
tnkm.pred = data.frame(tnkm = seq(from=0.74, to= 100, by= 0.1))
tnkm.pred$logit.percent.disturbed = -2.629818536 + -0.019183 *tnkm.pred$tnkm
tnkm.pred$logit.percent.disturbed.upper = -2.629818536 + (-0.019183+0.0051545) *tnkm.pred$tnkm
tnkm.pred$logit.percent.disturbed.lower = -2.629818536 + (-0.019183-0.0051545) *tnkm.pred$tnkm
tnkm.pred$percent.disturbed = (plogis(tnkm.pred$logit.percent.disturbed)-0.001)*100
tnkm.pred$percent.disturbed.upper = (plogis(tnkm.pred$logit.percent.disturbed.upper)-0.001)*100
tnkm.pred$percent.disturbed.lower = (plogis(tnkm.pred$logit.percent.disturbed.lower)-0.001)*100
plot(percent.disturbed~tnkm, data = tnkm.pred, xlab = “Distance to town (km)”, ylab=”Percentage of area disturbed”)
#start of plot
p = ggplot(tnkm.pred, aes(tnkm,percent.disturbed))
p + geom_line(colour=I(“blue”),size=1.1)+
geom_ribbon(aes(ymin=percent.disturbed.lower, ymax=percent.disturbed.upper),
fill=I(“blue”),alpha=0.2)+
theme_bw() +
labs(x=”Distance from town (kms)”,y=”Area disturbed (%)”)
ggsave(“dp figure 7 – distance from town.png”, width=4, height=3, dpi=300)
# end of plot

###########################################################################################################
# an afterthought: it would be good to combine pastoral and cons_estate as they both refer to tenure.
# I’ve resolved a few conflicts and created a new variable called ‘tenure’; the AIC for this model is slightly better,

dptenure = lme(ldistp ~ factor(tenure) + factor(esa) + factor(veg_form) + factor(greenstone) +
lntyrs + tnkm, random = ~1|SA_number,method = “REML”, data=dp)
summary(dptenure)
summary(dp$tenure)

# Maybe this would be better but I’d have to do all the predictions over again in excel and that’s too hard
# but for now I’ll have a sneak peak at the result of it:
# Class-A is the default,
# former leasehold has more disturbance
# gazetted cons has less disturbance
# not-cons has more disturbance but less than former leasehold
# pastoral has more disturbance but less than not-cons
# unofficial cons has more than pastoral

# so in order of most to least disturbance:
# 1. former leasehold
# 2. not-cons-estate
# 3. unofficial cons
# 4. pastoral
# 5. Class-A
# 6. gazetted conservation

# It’s wierd, you’d expect pastoral and former-pastoral to match up a bit more.
# ..and perhaps unofficial cons to match up more with gazetted cons.
# it looks as if the process of converting a pastoral lease to a conservation-ex-lease is the most disturbing thing!
# doesn’t really make sense, unless it’s all focued on Credo and they cleared a lot more area?
# I think i need to actually just combine the pastoral and ex-pastoral categories and that might increase the sample size and solve the issue.

Posted in blog | Leave a comment

A reflection

Two quotes from the recent newsletter of DBytes, a Newletter of the Environmental Decisions Group, seem to highlight a manifestation of Dr. Robert Sapolsky’s claim that the one trait that best defines and distinguishes us humans from other species is the ability to hold two contradictory ideas in our head, and yet continue on in the face of it:

“Many species are slipping away before we can even describe them. This IUCN Red List update shows that the scale of the global extinction crisis may be even greater than we thought.”   – Inger Andersen, Director General, IUCN press release

meanwhile, back in Australia…

“By 2019 Federal Environment Department spending will have declined by 38 per cent on 2013 levels while total government spending will have increased by 21.7 per cent. National spending on biodiversity has fallen to the lowest levels in more than a decade – with conservation getting just 5 cents in every 100 dollars spent by the government.”
– Matthew Rose, ACF press release

In the meantime, this recent flier from the Australian Greenspace Alliance reminds us that sustaining the natural environment is not just a ‘giving thing’: it helps to provide us with the things that make us well and happy:

a3-health-infographicv2-landscape

Posted in blog | Leave a comment

Presentation to the Environmental Protection Authority

I’ll be presenting the findings of my PhD research to the Western Australian Environmental Protection Authority (EPA) in early December 2016. The EPA has statutory obligations under the Environmental Protection Act 1986 to conduct environmental impact assessments, initiate measures to protect the environment from environmental harm and pollution and to provide advice to the Minister on environmental matters generally.

This will be an excellent opportunity to disseminate my research and influence policy and decisions relating to conservation of Western Australia’s considerable biological assets.

Keren Raiter_PhD flier.png

Posted in blog | Leave a comment

PHinisheD

I’m pleased to announce that the University of Western Australia has accepted my thesis as satisfying the requirements for the degree of Doctor of Philosophy. You can now call me Dr. 🙂

PhD front page

My thesis is now available online at http://research-repository.uwa.edu.au/files/10146368/THESIS_DOCTOR_OF  _PHILOSOPHY_RAITER_Keren_Gila_2016.pdf

Acknowledgements
בס”ד
One of the most incredible things about doing this PhD has been experiencing the sheer breadth and depth of generosity, support and interest that I have been so fortunate to receive along the way, and which have made this journey not only possible, but profoundly enriching. It is my hope that I have done all of it justice with this thesis and its associated outcomes and will continue to do so with the knowledge and experience that I now carry.
Firstly, I had the honour of receiving unwavering support from a diverse and distinguished team of supervisors, who were skilled at encouraging me to think broadly and explore the frontiers of what conservation needs, as well as always being present with my research even while they were stretched across the continent, and often the globe. Richard, I learnt so much from your wise advice and pragmatic approach to the unfolding of the research mystery. I also really appreciated the independence you gave me to follow my ‘researcher nose’, and your support to attend a number of very useful conferences. Suzanne, you have been a real role model for me, and your knowledge of the Great Western Woodlands and its stakeholders, as well as your thoughtful and encouraging feedback on my work have really benefited me. Hugh, your ability to turn a problem (“most of the impacts are too subtle to be visible”) into a conceptual framework that set the stage for my whole thesis was a huge gift for me, as was your encouragement to do conservation powerfully (“try to channel Attila the Hun or Napoleon on your rewrite”). Leonie, you were a later addition to my supervisory team and I really appreciated your enthusiasm and support for the predator investigation, as well as your help in developing the methods.
I feel extremely privileged to have conducted my field work in the largest and most intact
remaining temperate woodland on earth, work which would not have been possible without the funding, support, and assistance of many organisations and people. I am extremely grateful to the Wilderness Society for providing the funds for the field work and some of the motion-sensor cameras and to Gondwana Link for inspiring this work in the first place and then helping me to finish it with completion funding. The Great Western Woodlands Collaboration which consisted of Gondwana Link, The Wilderness Society, Pew Trusts and The Nature Conservancy also provided invaluable inspiration and information. I thank the University of Western Australia for bestowing me with the Robert and Maude Gledden Research Scholarship and providing me with its research infrastructure. I thank the School of Plant Biology for providing funds and support that assisted me to conduct my research, and the Centre of Excellence for Environmental Decisions for providing me with a top-up scholarship. I thank CSIRO for providing me with a studentship and a spacious and quiet desk with a great garden view that was the setting for many productive work days.
I feel fortunate to have been part of the Ecosystem Restoration and Intervention Ecology
research group at UWA, led by Richard Hobbs. Our weekly meetings, social gatherings and office camaraderie have been educational, supportive and fun. Thanks to Mandy Trueman, Cristina Ramalho, Bridget Johnson, Dawn Dickinson, Jodi Price, Heather Gordon, Rebecca Campbell, Tim Morald, Maggie Triska, Michael Wysong, Joanna Burgar, Rachel Standish, Melinda Moir, Mike Craig, Mike Perring, Jelena May, Todd Erickson, Erika Roper, Hillary Harrop, Christine Allen and others for being part of this. Thanks also to UWA academic support staff who have played a role in my academic development, especially Joanne Edmonston. Krystyna Haq, and Michael Azariadis. I have also been fortunate to be affiliated with the Land and Water team led by Suzanne Prober and others at CSIRO in Floreat, and thank Nat Raisbeck-Brown, Carl Gosper, Georg Wiehl, Garry Ogston, Anna Simonsen, Emma Woodward, Craig MacFarlane and others for your friendly inclusion, assistance and advice. In addition, I have been privileged to be part of a wider network of researchers in the Environmental Decisions Group, led by Hugh Possingham at University of Queensland, and thank this network for its insights and support.
Cliffs Natural Resources provided immense support and assistance on field trips. Thank you to Johnny Prefumo, Nicole Harry, Jeremy Shepherdson, Rob Howard, Kylie Wilkinson, Chris Dart, Neil Smith, Lorna McDonald and Belinda Madigan of the Cliffs Environment team who were incredibly helpful and responsive in this regard. The Department of Parks and Wildlife have also supported the fieldwork with approvals and licences, pertinent information and regular fire safely check-ins; thanks to Keith Morris, Julie Futter, Vanessa Jackson, Jennifer Jackson, David Algar, Sarah Comer and others. Thank you also to doggers Gordon Anderson and Stuart McEwan who were both very helpful, taught me to identify predator prints and scats and provided valuable insights into predator behaviour. Thanks also to Gorgeanna Story for doing the scat analysis.
Heartfelt appreciation also goes to Sue and Rolf Meeking who provided unique insights, local advice and wonderful hospitality both on their farm and in the bush, as well as logistical assistance which came during a difficult time and helped me to stay on track with my fieldwork. Other local contacts Coral Carter of Kalgoorlie, Rev. Dr Anna Killigrew and Rev. Peter Harrison of Koora Retreat, and John and Bernadette Cashmore also provided much appreciated hospitality along the way.
Discussions with Keith Bradby, Peter Price, Amanda Keesing, Wayne O’Sullivan, Megan Evans, Judith Harvey, Shapelle McNee, Peter-Jon Waddell, James O’Connor, Barry Traill, Mark Gardener, Charles Roche, Liz Fox and others throughout the course of this research helped me to better understand the conservation situation and needs of the GWW and direct my research focus. Workshops organised by Megan Evans and Kerrie Wilson of University of Queensland and Liz Fox and Cheryl Gole of Birdlife Australia were also immensely useful in this regard.
Much of the work presented in this thesis was performed with the assistance of a suite of volunteers. I wish to thank my intrepid field volunteers, Fiona Westcott, Kieran Golby, Stewart Bayford, Neal Birch, Ophir Levin, Bridget Johnson, Rebecca Campbell, Joanna Burgar and Michael Wysong for their eagerness to leave their creature comforts behind and join me on these remote trips. I also thank them for their interest in the research and for their many questions and ideas that helped me to develop the investigations; for their humour and good company that made these trips unforgettable, and for their hard work and dependability that ensured successful completion of the work and safe homeward returns. I also wish to thank the team of volunteers who helped me to digitize the disturbance footprint of vast swathes of the Great Western Woodlands over many hours: Ophir Levin, Julia Waite, Brad Desmond, and Rachel Omodei.
Thank you to other friends and colleagues Jaya Penelope, Krishna Rose, Ronen Steingold, Kiran Kigs, Kieran Golby, Naomi and Yonatan Li’el, Tegan Rourke, Cherie Carlo, Cristina Ramalho, Jocelyn Peyret, Asael Greenfeld, Ayala Ben Yosef, Judith Harvey, James O’Connor, Alison Hurst Emma Jack, Terry Farrell and others who have been so supportive and understanding and contributed to me having a wholesome life outside of my PhD. Thank you also to the other members of the Tealeaf Troubadours – Jaya Penelope, Jesse Williamson and Alex Hey and the poetry community, especially Jackson and Elio Novello for tolerating my busyness, engaging with my science poems and providing a creative outlet that has been an important part of the journey. Thank you to professionals Laura Harvey, Maria Arora and Hala Bitdorf who helped and taught me along the way.
Thanks to Tom Brooks, Peter Muirden, and Tim Sparks at Department of Water for their understanding and flexibility regarding work arrangements to facilitate field work and timely completion of my PhD. Thanks also to friends and colleagues at the Department; Gill White, Frances Miller, Renee Dixon, Shafiq Alam, Sue Tillman, Georgina Evans, Jaci Moore and others for their interest and encouragement.
Finally, I heartily thank my family for their support. Thank you to my amazing parents, Raul and Lesley, who taught me that I can do anything I wanted and then supported me morally and practically to do so, especially with nourishing food and funds when the going got tough. Thanks to my sister, Perla, for your understanding, the use of your car, and those wonderful Shabbat dinners that were the best antidote to a long hard week. Thank you to my lively little nieces, Izabell, Katana and Scarlett, who provided embodied reasons to work towards a more ecologically harmonious future, and persuaded me to play even when I hadn’t finished that manuscript. Thanks also to my family around the world for their love, encouragement, and understanding, especially the Abuahrons, the Cheniks, and the Glazners. Mbapani Ngitoria was a great support in my process of embarking on this PhD journey and choosing the project that I would pursue. Sunny Blundell-Wignall held my hand through the last part of this journey, and gave me steadfast support as well as helping me to balance work with rest and play. I thank him and also his family for putting up with my absences and being so understanding.

 

Doctor.JPG

Posted in blog | 1 Comment

PhD completion seminar – March 14 2016

It is with great pleasure that I invite all interested people to attend my PhD completion seminar on Monday 14th March. The seminar will be a one-hour overview of my PhD research on the enigmatic ecological impacts of mining and linear infrastructure in the Great Western Woodlands of southwestern Australia.

When: 4pm Monday 14th March 2016

Where: Botany Seminar Room, Botany Building, The University of Western Australia

Why: To share the insights and results emerging from a mammoth research project that I have undertaken over the last four years.

Keren Raiter PhD completion seminar

Posted in blog | 1 Comment

Art and inspiration, Bungalbin

I was recently invited to contribute to the Wilderness Society’s exhibition of art inspired by the Great Western Woodlands, as part of an event held to celebrate the incredible but threatened Helena-Aurora Range, and gather support for its protection. To learn more about the event, the need for protection of Helena-Aurora Range, and ways of supporting its protection, click here.

The images and text that I exhibited and performed follow.

many-limbed

Many-limbed. Photograph © Keren Gila Raiter

Sunrise over Helena-Aurora. Photograph by Fiona Westcott (reproduced with permission)

Sunrise over Helena-Aurora. Photograph © Fiona Westcott (reproduced with permission)

Dianella revoluta – blueberry lilly. Photograph © Keren Gila Raiter

Moonrise almost as old as the moon fresh as the last burn

Moonrise

almost as old as the moon
fresh as the last burn

2. Verticordia chrysantha

Verticordia chrysantha

This is not the middle of nowhere
it’s the centre of everywhere
the sweet space between wet forests and dry deserts
where there’s more eucalypts than there’s elements in the periodic table;
more flowering plants than in the UK
where banded iron never goes out of fashion
with water in rocky cracks and rare views over subtle topography
what’s more, it’s my home

3. Bungalbin

Bungalbin

by the time we reached Bungalbin
we had forgotten what a hill looks like
and a range seemed impossible in this flat expansiveness
but the earth reaching skyward was unmistakable

We camped in Helena and Aurora’s wide embrace
long ironstone arms stretched out around us
striped with geology

4. Night creature

Night lacewing (Myrmeleontidae family)

the creatures of the night
remind me of the mystery
of life
of ecosystems
of the things that are hidden from our view
but that are nevertheless
essential parts of our existence

Abandoned mine

Abandoned mine

after the minerals have been traded
profits spent
workers retrenched, or retired
that water will still be a strange shade of green

6. red legged- arachnid

Red-legged arachnid

I walked a thousand kilometres
till my legs were red and hairy
I lost and found myself
in between these leaves and branches
and I won’t forget

It's us who decide 1 Slide2 Slide3 Slide4

All photographs and text © Keren Gila Raiter except where noted

Posted in blog | 6 Comments