Chapter 21 Generalized linear model
21.1 Learning goals
- Logistic regression.
- Logit transform.
- Fitting a logistic regression in R.
- Visualizing and interpreting model predictions.
- Simulating data from a logistic regression.
- Assessing model fit.
- Testing hypotheses.
- Reporting results.
- Mixed effects logistic regression.
21.2 Load packages and set plotting theme
library("knitr") # for knitting RMarkdown
library("titanic") # titanic dataset
library("kableExtra") # for making nice tables
library("janitor") # for cleaning column names
library("broom") # for tidying up linear models
library("lme4") # for linear mixed effects models
library("boot") # for bootstrapping (also has an inverse logit function)
library("ggeffects") # for showing effects in linear, generalized linear, and other models
library("afex") # for significance testing of mixed effects models
library("emmeans") # for the joint_tests() function
library("tidyverse") # for wrangling, plotting, etc.
theme_set(theme_classic() + #set the theme
theme(text = element_text(size = 20))) #set the default text size
$set(comment = "",
opts_chunkfig.show = "hold")
21.3 Load data set
= titanic_train %>%
df.titanic clean_names() %>%
mutate(sex = as.factor(sex))
Let’s take a quick look at the data:
%>%
df.titanic glimpse()
Rows: 891
Columns: 12
$ passenger_id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ survived <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, …
$ pclass <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3, …
$ name <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley (F…
$ sex <fct> male, female, female, female, male, male, male, male, fem…
$ age <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 14,…
$ sib_sp <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, …
$ parch <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0, …
$ ticket <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", "3…
$ fare <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625…
$ cabin <chr> "", "C85", "", "C123", "", "", "E46", "", "", "", "G6", "…
$ embarked <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", "S…
# Table of the first 10 entries
%>%
df.titanic head(10) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped",
full_width = F)
| passenger_id | survived | pclass | name | sex | age | sib_sp | parch | ticket | fare | cabin | embarked |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22 | 1 | 0 | A/5 21171 | 7.25 | S | |
| 2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Thayer) | female | 38 | 1 | 0 | PC 17599 | 71.28 | C85 | C |
| 3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26 | 0 | 0 | STON/O2. 3101282 | 7.92 | S | |
| 4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35 | 1 | 0 | 113803 | 53.10 | C123 | S |
| 5 | 0 | 3 | Allen, Mr. William Henry | male | 35 | 0 | 0 | 373450 | 8.05 | S | |
| 6 | 0 | 3 | Moran, Mr. James | male | NA | 0 | 0 | 330877 | 8.46 | Q | |
| 7 | 0 | 1 | McCarthy, Mr. Timothy J | male | 54 | 0 | 0 | 17463 | 51.86 | E46 | S |
| 8 | 0 | 3 | Palsson, Master. Gosta Leonard | male | 2 | 3 | 1 | 349909 | 21.08 | S | |
| 9 | 1 | 3 | Johnson, Mrs. Oscar W (Elisabeth Vilhelmina Berg) | female | 27 | 0 | 2 | 347742 | 11.13 | S | |
| 10 | 1 | 2 | Nasser, Mrs. Nicholas (Adele Achem) | female | 14 | 1 | 0 | 237736 | 30.07 | C |
21.4 Logistic regression
Let’s see if we can predict whether or not a passenger survived based on the price of their ticket.
Let’s run a simple regression first:
# fit a linear model
= lm(formula = survived ~ 1 + fare,
fit.lm data = df.titanic)
# summarize the results
%>%
fit.lm summary()
Call:
lm(formula = survived ~ 1 + fare, data = df.titanic)
Residuals:
Min 1Q Median 3Q Max
-0.9653 -0.3391 -0.3222 0.6044 0.6973
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3026994 0.0187849 16.114 < 2e-16 ***
fare 0.0025195 0.0003174 7.939 6.12e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4705 on 889 degrees of freedom
Multiple R-squared: 0.06621, Adjusted R-squared: 0.06516
F-statistic: 63.03 on 1 and 889 DF, p-value: 6.12e-15
Look’s like fare
is a significant predictor of whether or not a person survived. Let’s visualize the model’s predictions:
ggplot(data = df.titanic,
mapping = aes(x = fare,
y = survived)) +
geom_smooth(method = "lm") +
geom_point() +
labs(y = "survived")
This doesn’t look good! The model predicts intermediate values of survived
(which doesn’t make sense given that a person either survived or didn’t survive). Furthermore, the model predicts values greater than 1 for fares greather than ~ 300.
Let’s run a logistic regression instead.
# fit a logistic regression
= glm(formula = survived ~ 1 + fare,
fit.glm family = "binomial",
data = df.titanic)
%>%
fit.glm summary()
Call:
glm(formula = survived ~ 1 + fare, family = "binomial", data = df.titanic)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4906 -0.8878 -0.8531 1.3429 1.5942
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.941330 0.095129 -9.895 < 2e-16 ***
fare 0.015197 0.002232 6.810 9.79e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1186.7 on 890 degrees of freedom
Residual deviance: 1117.6 on 889 degrees of freedom
AIC: 1121.6
Number of Fisher Scoring iterations: 4
And let’s visualize the predictions of the logistic regression:
ggplot(data = df.titanic,
mapping = aes(x = fare,
y = survived)) +
geom_smooth(method = "glm",
method.args = list(family = "binomial")) +
geom_point() +
labs(y = "p(survived)")
Much better! Note that we’ve changed the interpretation of our dependent variable. We are now predicting the probability that a person survived based on their fare. The model now only predicts values between 0 and 1. To achieve this, we apply a logit transform to the outcome variable like so:
\[ \ln(\frac{\pi_i}{1-\pi_i}) = b_0 + b_1 \cdot X_i + e_i \] where \(\pi_i\) is the probability of passenger \(i\) having survived. Importantly, this affects our interpretation of the model parameters. They are now defined in log-odds, and can apply an inverse logit transformation to turn this back into a probability:
With
\[ \pi = P(Y = 1) \] and the logit transformation
\[ \ln(\frac{\pi}{1-\pi}) = V, \] where \(V\) is just a placeholder for our linear model formula, we can go back to \(\pi\) through the inverse logit transformation like so:
\[
\pi = \frac{e^V}{1 + e^V}
\]
In R, we can use log(x)
to calculate the natural logarithm \(\ln(x)\), and exp(x)
to calculate e^x
.
21.4.1 Interpreting the parameters
%>%
fit.glm summary()
Call:
glm(formula = survived ~ 1 + fare, family = "binomial", data = df.titanic)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4906 -0.8878 -0.8531 1.3429 1.5942
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.941330 0.095129 -9.895 < 2e-16 ***
fare 0.015197 0.002232 6.810 9.79e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1186.7 on 890 degrees of freedom
Residual deviance: 1117.6 on 889 degrees of freedom
AIC: 1121.6
Number of Fisher Scoring iterations: 4
The estimate for the intercept and fare are in log-odds.
Let’s take a look at the linear model’s predictions in log-odds space.
= fit.glm %>%
df.plot augment() %>%
clean_names()
ggplot(data = df.plot,
mapping = aes(x = fare,
y = fitted)) +
geom_line()
Nice, looks like a good old linear model. But what’s the y-axis here? It’s in log-odds (buh!). This is difficult to interpret. Let’s transform the y-axis back to probabilities to make it easier to see what’s going on.
ggplot(data = df.plot,
mapping = aes(x = fare,
y = inv.logit(fitted))) +
geom_line()
Great! Now the y-axis is back in probability space. We used the inverse logit function inv.logit()
to transfer to log-odds back into probabilities.
Let’s check what the intercept of our model is in probability space:
%>%
fit.glm pluck(coefficients, 1) %>%
inv.logit()
[1] 0.2806318
Here, we see that the intercept is \(p = 0.28\). That is, the predicted chance of survival for someone who didn’t pay any fare at all is 28% according to the model. Interpreting the slope is a little more tricky. Let’s look at a situation first where we have a binary predictor.
<IP_ADDRESS> Logit transform
Here is a visualization of what the odds and log odds transformation look like.
# going from probabilities to odds (ranges from 0 to +infinity)
ggplot(data = tibble(x = seq(0, 1, 0.1)),
mapping = aes(x = x)) +
stat_function(fun = ~ ./(1 - .),
size = 1)
# going from probabilities to log odds (ranges from -infinity to +infinity)
ggplot(data = tibble(x = seq(0, 1, 0.1)),
mapping = aes(x = x)) +
stat_function(fun = ~ log(./(1 - .)),
size = 1)
21.4.2 Binary predictor
Let’s see whether the probability of survival differed between male and female passengers.
= glm(formula = survived ~ 1 + sex,
fit.glm2 family = "binomial",
data = df.titanic)
%>%
fit.glm2 summary()
Call:
glm(formula = survived ~ 1 + sex, family = "binomial", data = df.titanic)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6462 -0.6471 -0.6471 0.7725 1.8256
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.0566 0.1290 8.191 2.58e-16 ***
sexmale -2.5137 0.1672 -15.036 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1186.7 on 890 degrees of freedom
Residual deviance: 917.8 on 889 degrees of freedom
AIC: 921.8
Number of Fisher Scoring iterations: 4
It looks like it did! Let’s visualize:
%>%
df.titanic mutate(survived = factor(survived, labels = c("died", "survived"))) %>%
ggplot(data = .,
mapping = aes(x = sex,
fill = survived)) +
geom_bar(position = "fill",
color = "black") +
scale_fill_brewer(palette = "Set1") +
labs(x = "",
fill = "",
y = "probability")
And let’s interpret the parameters by applying the inverse logit transform. To get the prediction for female passengers we get
\[ \widehat{\ln(\frac{\pi_i}{1-\pi_i})} = b_0 + b_1 \cdot \text{sex}_i = b_0 + b_1 \cdot 0 = b_0 \] since we dummy coded the predictor and female is our reference category. To get the predicted probability of survival for women we do the logit transform:
\[ \pi = \frac{e^{b_0}}{1 + e^{b_0}} \] The predicted probability is:
%>%
fit.glm2 pluck(coefficients, 1) %>%
inv.logit()
[1] 0.7420382
To get the prediction for male passengers we have:
\[ \widehat{\ln(\frac{\pi_i}{1-\pi_i})} = b_0 + b_1 \cdot \text{sex}_i = b_0 + b_1 \cdot 1 = b_0 + b_1 \] Applying the logit transform like so
\[ \pi = \frac{e^{b_0 + b_1}}{1 + e^{b_0 + b_1}} \]
The predicted probability of male passengers surviving is:
%>%
fit.glm2 pluck(coefficients) %>%
sum() %>%
inv.logit()
[1] 0.1889081
Here is the same information in a table:
%>%
df.titanic count(sex, survived) %>%
mutate(p = n / sum(n)) %>%
group_by(sex) %>%
mutate(`p(survived|sex)` = p / sum(p))
# A tibble: 4 × 5
# Groups: sex [2]
sex survived n p `p(survived|sex)`
<fct> <int> <int> <dbl> <dbl>
1 female 0 81 0.0909 0.258
2 female 1 233 0.262 0.742
3 male 0 468 0.525 0.811
4 male 1 109 0.122 0.189
21.4.3 Continuous predictor
To interpret the predictions when a continuous predictor is involved, it’s easiest to consider a few concrete cases. Here, I use the augment()
function from the “broom” package to get the model’s predictions for some values of interest:
%>%
fit.glm augment(newdata = tibble(fare = c(0, 10, 50, 100, 500))) %>%
clean_names() %>%
select(fare, prediction = fitted) %>%
mutate(`p(survival)` = inv.logit(prediction))
# A tibble: 5 × 3
fare prediction `p(survival)`
<dbl> <dbl> <dbl>
1 0 -0.941 0.281
2 10 -0.789 0.312
3 50 -0.181 0.455
4 100 0.578 0.641
5 500 6.66 0.999
21.4.4 Several predictors
Let’s fit a logistic regression that predicts the probability of survival based both on the passenger’s sex and what fare they paid (allowing for an interaction of the two predictors):
= glm(formula = survived ~ 1 + sex * fare,
fit.glm3 family = "binomial",
data = df.titanic)
%>%
fit.glm3 summary()
Call:
glm(formula = survived ~ 1 + sex * fare, family = "binomial",
data = df.titanic)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6280 -0.6279 -0.5991 0.8172 1.9288
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.408428 0.189999 2.150 0.031584 *
sexmale -2.099345 0.230291 -9.116 < 2e-16 ***
fare 0.019878 0.005372 3.701 0.000215 ***
sexmale:fare -0.011617 0.005934 -1.958 0.050252 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1186.66 on 890 degrees of freedom
Residual deviance: 879.85 on 887 degrees of freedom
AIC: 887.85
Number of Fisher Scoring iterations: 5
Make sure not to interpret the significance test on the coefficients here as main effects. Based on this summary table, you cannot say whether there is a significant difference between male vs. female passenger in their probability of survival. What coefficient for sexmale
captures is whether there is a significant difference between male and female passengers who paid a fare of 0. That is, it’s the predicted difference between the reference category (female) and the other category (male) when all other predictors are 0.
Let’s visualize the model predictions:
ggplot(data = df.titanic,
mapping = aes(x = fare,
y = survived,
color = sex)) +
geom_point(alpha = 0.1, size = 2) +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
alpha = 0.2,
aes(fill = sex)) +
scale_color_brewer(palette = "Set1")
Just for kicks, to get a better sense for what the interaction looks like, here is the visualization in log-odds space:
%>%
fit.glm3 augment() %>%
clean_names() %>%
ggplot(data = .,
mapping = aes(x = fare,
color = sex,
y = fitted)) +
geom_line() +
scale_color_brewer(palette = "Set1")
Let’s see how large the difference between genders is once we take into account how much each person paid for the fair:
ggpredict(fit.glm3,
terms = c("sex"))
# Predicted probabilities of survived
sex | Predicted | 95% CI
---------------------------------
female | 0.74 | [0.69, 0.79]
male | 0.19 | [0.16, 0.23]
Adjusted for:
* fare = 32.20
We notice that there is one outlier who was male and paid a $500 fare (or maybe this is a mistake in the data entry?!). Let’s remove this outlier and see what happens:
= glm(formula = survived ~ 1 + sex * fare,
fit.glm3_no_outlier family = "binomial",
data = df.titanic %>%
filter(fare < 500))
%>%
fit.glm3_no_outlier summary()
Call:
glm(formula = survived ~ 1 + sex * fare, family = "binomial",
data = df.titanic %>% filter(fare < 500))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6280 -0.6288 -0.6018 0.8172 1.9226
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.408436 0.190019 2.149 0.031598 *
sexmale -2.085344 0.232260 -8.979 < 2e-16 ***
fare 0.019878 0.005373 3.699 0.000216 ***
sexmale:fare -0.012178 0.006066 -2.008 0.044688 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1180.89 on 887 degrees of freedom
Residual deviance: 879.51 on 884 degrees of freedom
AIC: 887.51
Number of Fisher Scoring iterations: 5
%>%
df.titanic filter(fare < 500) %>%
mutate(sex = as.factor(sex)) %>%
ggplot(data = .,
mapping = aes(x = fare,
y = survived,
color = sex)) +
geom_point(alpha = 0.1, size = 2) +
stat_smooth(method = "glm",
method.args = list(family = "binomial"),
alpha = 0.2,
fullrange = T,
aes(fill = sex)) +
scale_color_brewer(palette = "Set1") +
scale_x_continuous(limits = c(0, 500))
There is still a clear difference between female and male passengers, but the prediction for male passengers has changed a bit. Let’s look at a concrete example:
# with the outlier:
# predicted probability of survival for a male passenger who paid $200 for their fare
inv.logit(fit.glm3$coefficients[1] + fit.glm3$coefficients[2] +
$coefficients[3] * 200 + fit.glm3$coefficients[4] * 200) fit.glm3
(Intercept)
0.4903402
# without the outlier:
# predicted probability of survival for a male passenger who paid $200 for their fare
inv.logit(fit.glm3_no_outlier$coefficients[1] +
$coefficients[2] +
fit.glm3_no_outlier$coefficients[3] * 200 +
fit.glm3_no_outlier$coefficients[4] * 200) fit.glm3_no_outlier
(Intercept)
0.4658284
With the outlier removed, the predicted probability of survival for a male passenger who paid $200 decreases from 49% to 47%.
21.4.5 Using the “ggeffects” package
The “ggeffects” package helps with the interpretation of the results. It applies the inverse logit transform for us, and shows the predictions for a range of cases.
# show effects
ggeffect(model = fit.glm,
terms = "fare [0, 100, 200, 300, 400, 500]")
# Predicted probabilities of survived
fare | Predicted | 95% CI
-------------------------------
0 | 0.28 | [0.24, 0.32]
100 | 0.64 | [0.56, 0.72]
200 | 0.89 | [0.79, 0.95]
300 | 0.97 | [0.92, 0.99]
400 | 0.99 | [0.97, 1.00]
500 | 1.00 | [0.99, 1.00]
I’ve used the []
construction to specify for what values of the predictor fare
, I’d like get the predicted values. Here, the prediction is marginalized across both women and men.
We can also get a plot of the model predictions like so:
ggeffect(model = fit.glm,
terms = "fare [0, 100, 200, 300, 400, 500]") %>%
plot()
And, we can also get the predicted probability of survival for sex marginalized over the fare, using the model which included both sex and fare, as well as its interaction as predictors.
ggeffect(model = fit.glm3,
terms = "sex")
# Predicted probabilities of survived
sex | Predicted | 95% CI
---------------------------------
female | 0.74 | [0.69, 0.79]
male | 0.19 | [0.16, 0.23]
Finally, we can ask for predictions for specific combinations of our predictor variables, by using the ggpredict()
function.
ggpredict(model = fit.glm3,
terms = c("sex", "fare [200]"))
# Predicted probabilities of survived
sex | Predicted | 95% CI
---------------------------------
female | 0.99 | [0.93, 1.00]
male | 0.49 | [0.29, 0.70]
The example above, shows the predicted probability of survival for male vs. female passengers, assuming that they paid 200 for the fare.
21.5 Simulate a logistic regression
As always, to better understand a statistical modeling procedure, it’s helpful to simulate data from the assumed data-generating process, fit the model, and see whether we can reconstruct the parameters.
# make example reproducible
set.seed(1)
# set parameters
= 1000
sample_size = 0
b0 = 1
b1 # b1 = 8
# generate data
= tibble(x = rnorm(n = sample_size),
df.data y = b0 + b1 * x,
p = inv.logit(y)) %>%
mutate(response = rbinom(n(), size = 1, p = p))
# fit model
= glm(formula = response ~ 1 + x,
fit family = "binomial",
data = df.data)
# model summary
%>%
fit summary()
Call:
glm(formula = response ~ 1 + x, family = "binomial", data = df.data)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1137 -1.0118 -0.4591 1.0287 2.2591
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.06214 0.06918 -0.898 0.369
x 0.92905 0.07937 11.705 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1385.4 on 999 degrees of freedom
Residual deviance: 1209.6 on 998 degrees of freedom
AIC: 1213.6
Number of Fisher Scoring iterations: 3
%>%
df.data head(10) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped",
full_width = F)
| x | y | p | response |
|---|---|---|---|
| -0.63 | -0.63 | 0.35 | 1 |
| 0.18 | 0.18 | 0.55 | 0 |
| -0.84 | -0.84 | 0.30 | 1 |
| 1.60 | 1.60 | 0.83 | 1 |
| 0.33 | 0.33 | 0.58 | 1 |
| -0.82 | -0.82 | 0.31 | 0 |
| 0.49 | 0.49 | 0.62 | 1 |
| 0.74 | 0.74 | 0.68 | 1 |
| 0.58 | 0.58 | 0.64 | 1 |
| -0.31 | -0.31 | 0.42 | 0 |
Nice! The inferred estimates are very close to the parameter values we used to simulate the data.
Let’s visualize the result:
ggplot(data = df.data,
mapping = aes(x = x,
y = response)) +
geom_smooth(method = "glm",
method.args = list(family = "binomial")) +
geom_point(alpha = 0.1) +
labs(y = "p(response)")
21.5.1 Calculate the model’s likelihood
To calculate the likelihood of the data for a given logistic model, we look at the actual response, and the probability of the predicted response, and then determine the likelihood of the observation assuming a Bernoulli process. To get the overall likelihood of the data, we then multiply the likelihood of each data point (or take the logs first and then the sum to get the log-likelihood).
This table illustrate the steps involved:
%>%
fit augment() %>%
clean_names() %>%
mutate(p = inv.logit(fitted)) %>%
select(response, p) %>%
mutate(p_response = ifelse(response == 1, p, 1-p),
log_p = log(p_response)) %>%
rename(`p(Y = 1)` = p, `p(Y = response)` = p_response,
`log(p(Y = response))` = log_p)
# A tibble: 1,000 × 4
response `p(Y = 1)` `p(Y = response)` `log(p(Y = response))`
<int> <dbl> <dbl> <dbl>
1 1 0.344 0.344 -1.07
2 0 0.527 0.473 -0.749
3 1 0.302 0.302 -1.20
4 1 0.805 0.805 -0.217
5 1 0.561 0.561 -0.579
6 0 0.305 0.695 -0.364
7 1 0.596 0.596 -0.517
8 1 0.651 0.651 -0.429
9 1 0.616 0.616 -0.484
10 0 0.414 0.586 -0.535
# … with 990 more rows
Let’s calculate the log-likelihood by hand:
%>%
fit augment() %>%
clean_names() %>%
mutate(p = inv.logit(fitted),
log_likelihood = response * log(p) + (1 - response) * log(1 - p)) %>%
summarize(log_likelihood = sum(log_likelihood))
# A tibble: 1 × 1
log_likelihood
<dbl>
1 -605.
And compare it with the model summary
%>%
fit glance() %>%
select(logLik, AIC, BIC)
# A tibble: 1 × 3
logLik AIC BIC
<dbl> <dbl> <dbl>
1 -605. 1214. 1223.
We’re getting the same result – neat! Now we know how the likelihood of the data is calculated for a logistic regression model.
21.6 Testing hypotheses
To test hypotheses, we can use our good old model comparison approach:
# fit compact model
= glm(formula = survived ~ 1 + fare,
fit.compact family = "binomial",
data = df.titanic)
# fit augmented model
= glm(formula = survived ~ 1 + sex + fare,
fit.augmented family = "binomial",
data = df.titanic)
# likelihood ratio test
anova(fit.compact, fit.augmented, test = "LRT")
Analysis of Deviance Table
Model 1: survived ~ 1 + fare
Model 2: survived ~ 1 + sex + fare
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 889 1117.57
2 888 884.31 1 233.26 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note that in order to get a p-value out of this, we need to specify what statistical test we’d like to run. In this case, we use the likelihood ratio test (“LRT”).
We can also test for both effects of survived and fare in one go using the joint_tests()
function from the “emmeans” package like so:
glm(formula = survived ~ 1 + sex + fare,
family = "binomial",
data = df.titanic) %>%
joint_tests()
model term df1 df2 F.ratio p.value
sex 1 Inf 201.881 <.0001
fare 1 Inf 23.869 <.0001
Notice that the F.ratio reported using joint_tests()
(201.881) is not quite the same as the deviance value that we get through the likelihood ratio test (233.26). You can read more about why these two can come apart here.
21.7 Logistic mixed effects model
Just like we can build linear mixed effects models using lmer()
instead of lm()
, we can also build a logistic mixed effects regression using glmer()
instead of glm()
.
Let’s read in some data:
# load bdf data set from nlme package
data(bdf, package = "nlme")
= bdf %>%
df.language clean_names() %>%
filter(repeatgr != 2) %>%
mutate(repeatgr = repeatgr %>%
as.character() %>%
as.numeric())
rm(bdf)
Fit the model, and print out the results:
= glmer(repeatgr ~ 1 + ses + minority + (1 | school_nr),
fit data = df.language,
family = "binomial")
%>%
fit summary()
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: repeatgr ~ 1 + ses + minority + (1 | school_nr)
Data: df.language
AIC BIC logLik deviance df.resid
1659.1 1682.1 -825.6 1651.1 2279
Scaled residuals:
Min 1Q Median 3Q Max
-0.9235 -0.4045 -0.3150 -0.2249 5.8372
Random effects:
Groups Name Variance Std.Dev.
school_nr (Intercept) 0.2489 0.4989
Number of obs: 2283, groups: school_nr, 131
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.506291 0.197566 -2.563 0.01039 *
ses -0.060086 0.007523 -7.986 1.39e-15 ***
minorityY 0.673612 0.238657 2.823 0.00476 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) ses
ses -0.898
minorityY -0.308 0.208
To visualize the results, we can use the ggeffects
package.
ggpredict(model = fit,
terms = c("ses [all]", "minority")) %>%
plot()
And for significance testing, we can use the the joint_tests()
function from the “emmeans” package
glmer(formula = repeatgr ~ 1 + ses + minority + (1 | school_nr),
data = df.language,
family = "binomial") %>%
joint_tests()
model term df1 df2 F.ratio p.value
ses 1 Inf 63.784 <.0001
minority 1 Inf 7.967 0.0048
The results show that there was both a significant effect of ses and of minority.
Note: This post here says a little more about the relationship of the F.ratio in the joint_tests()
function, and what a likelihood ratio test yields. In short, it’s roughly the same thing.
If you’d like to compute the likelihood ratio test, a convenient way of doing so is by using the mixed()
function from the “afex” package.
mixed(formula = repeatgr ~ 1 + ses + minority + (1 | school_nr),
family = "binomial",
data = df.language,
method = "LRT")
Contrasts set to contr.sum for the following variables: minority, school_nr
Numerical variables NOT centered on 0: ses
If in interactions, interpretation of lower order (e.g., main) effects difficult.
Mixed Model Anova Table (Type 3 tests, LRT-method)
Model: repeatgr ~ 1 + ses + minority + (1 | school_nr)
Data: df.language
Df full model: 4
Effect df Chisq p.value
1 ses 1 75.39 *** <.001
2 minority 1 7.53 ** .006
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
And we can compare that the model comparison approach gives us the same result:
= glmer(repeatgr ~ 1 + ses + minority + (1 | school_nr),
fit_a data = df.language,
family = "binomial")
# dropping ses as a predictor
= glmer(repeatgr ~ 1 + minority + (1 | school_nr),
fit_c data = df.language,
family = "binomial")
anova(fit_a, fit_c, test = "LRT")
Data: df.language
Models:
fit_c: repeatgr ~ 1 + minority + (1 | school_nr)
fit_a: repeatgr ~ 1 + ses + minority + (1 | school_nr)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_c 3 1732.5 1749.7 -863.27 1726.5
fit_a 4 1659.1 1682.1 -825.57 1651.1 75.395 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
21.9 Session info
Information about this R session including which version of R was used, and what packages were loaded.
sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
[5] readr_2.1.2 tidyr_1.2.0 tibble_3.1.7 ggplot2_3.3.6
[9] tidyverse_1.3.1 emmeans_1.7.3 afex_1.1-1 ggeffects_1.1.2
[13] boot_1.3-28 lme4_1.1-29 Matrix_1.4-1 broom_0.8.0
[17] janitor_2.1.0 kableExtra_1.3.4 titanic_0.1.0 knitr_1.39
loaded via a namespace (and not attached):
[1] TH.data_1.1-1 minqa_1.2.4 colorspace_2.0-3
[4] ellipsis_0.3.2 sjlabelled_1.2.0 estimability_1.3
[7] snakecase_0.11.0 fs_1.5.2 rstudioapi_0.13
[10] farver_2.1.0 fansi_1.0.3 mvtnorm_1.1-3
[13] lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18
[16] splines_4.1.2 effects_4.2-1 jsonlite_1.8.0
[19] nloptr_2.0.0 dbplyr_2.1.1 compiler_4.1.2
[22] httr_1.4.3 backports_1.4.1 assertthat_0.2.1
[25] fastmap_1.1.0 survey_4.1-1 cli_3.3.0
[28] htmltools_0.5.2 tools_4.1.2 lmerTest_3.1-3
[31] coda_0.19-4 gtable_0.3.0 glue_1.6.2
[34] reshape2_1.4.4 Rcpp_<IP_ADDRESS> carData_3.0-5
[37] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.4.1
[40] svglite_2.1.0 nlme_3.1-157 insight_0.17.0
[43] xfun_0.30 rvest_1.0.2 lifecycle_1.0.1
[46] MASS_7.3-57 zoo_1.8-10 scales_1.2.0
[49] hms_1.1.1 parallel_4.1.2 sandwich_3.0-1
[52] RColorBrewer_1.1-3 yaml_2.3.5 sass_0.4.1
[55] stringi_1.7.6 highr_0.9 rlang_1.0.2
[58] pkgconfig_2.0.3 systemfonts_1.0.4 evaluate_0.15
[61] lattice_0.20-45 labeling_0.4.2 tidyselect_1.1.2
[64] plyr_1.8.7 magrittr_2.0.3 bookdown_0.26
[67] R6_2.5.1 generics_0.1.2 multcomp_1.4-19
[70] DBI_1.1.2 pillar_1.7.0 haven_2.5.0
[73] withr_2.5.0 mgcv_1.8-40 survival_3.3-1
[76] abind_1.4-5 nnet_7.3-17 modelr_0.1.8
[79] crayon_1.5.1 car_3.0-13 utf8_1.2.2
[82] tzdb_0.3.0 rmarkdown_2.14 grid_4.1.2
[85] readxl_1.4.0 reprex_2.0.1 digest_0.6.29
[88] webshot_0.5.3 xtable_1.8-4 numDeriv_2016.8-1.1
[91] munsell_0.5.0 viridisLite_0.4.0 bslib_0.3.1
[94] mitools_2.4


Better than email for organizing pickup basketball, soccer, and ultimate.
{{#unless SGet 'sign-in'}} Add your info or sign in {{else}} Sign in or add your info {{/unless}} to {{action}}.
Add friend
Check your email for more info.
{{> alerts where="subscribe"}}Would you like to receive email when new {{type}} games are posted near {{placeName}}?
{{> subscribeButton}}You can use search to find locations and game types to which you want to subscribe.
{{#if loading}} loading your subscriptions... {{else}} {{#if hasSubs}}You have {{subsCount}} subscriptions.
This page will get fancier soon, enabling you to view and edit subscriptions.
{{else}} You have no subscriptions. {{/if}} {{/if}}
{{#each comments}}- {{> gameComment}}
{{/each}}


Predictive models of malaria vector larval habitat locations may provide a basis for understanding the spatial determinants of malaria transmission.
We used four landscape variables (topographic wetness index [TWI], soil type, land use-land cover, and distance to stream) and accumulated precipitation to model larval habitat locations in a region of western Kenya through two methods: logistic regression and random forest. Additionally, we used two separate data sets to account for variation in habitat locations across space and over time.
Larval habitats were more likely to be present in locations with a lower slope to contributing area ratio (i.e. TWI), closer to streams, with agricultural land use relative to nonagricultural land use, and in friable clay/sandy clay loam soil and firm, silty clay/clay soil relative to friable clay soil. The probability of larval habitat presence increased with increasing accumulated precipitation. The random forest models were more accurate than the logistic regression models, especially when accumulated precipitation was included to account for seasonal differences in precipitation. The most accurate models for the two data sets had area under the curve (AUC) values of 0.864 and 0.871, respectively. TWI, distance to the nearest stream, and precipitation had the greatest mean decrease in Gini impurity criteria in these models.
This study demonstrates the usefulness of random forest models for larval malaria vector habitat modeling. TWI and distance to the nearest stream were the two most important landscape variables in these models. Including accumulated precipitation in our models improved the accuracy of larval habitat location predictions by accounting for seasonal variation in the precipitation. Finally, the sampling strategy employed here for model parameterization could serve as a framework for creating predictive larval habitat models to assist in larval control efforts.
Malaria is one of the most significant infectious diseases affecting people in poverty, with an estimated 219 million cases of malaria worldwide in 2010 killing 660,000 people [ 1]. An estimated 1.44 billion people in South America, Africa, and Asia lived in areas with stable transmission of malaria caused by Plasmodium falciparum (Welch), yet the risk of P. falciparum transmission varies considerably across its range [ 2]. Even at fine scales, the spatial distribution of malaria is heterogeneous, differing among households within a community [ 3- 6]. Of course, socioeconomic and immunological differences contribute to the spatial heterogeneity of malaria [ 3, 4]. Additionally, landscape factors contribute to the spatial distribution of malaria [ 5, 7], which is likely an indirect relationship ultimately due, in part, to the influence of landscape factors on the locations of the aquatic habitats of the vector mosquito larvae. The spatial distribution of the larval habitats partially determines the spatial distribution of the adult malaria vectors in many landscapes [ 5, 8- 10]. Subsequently, the heterogeneous spatial distribution of malaria vectors among households coincides with the spatial distribution of malaria parasitemia in some landscapes [ 5]. Therefore, understanding the factors that determine the distribution of the larval habitats facilitates our understanding of the spatial determinants of malaria transmission.
The vast majority of deaths from malaria (91%) occur in Africa [ 1], where the primary vector mosquitoes are among the most efficient vectors of malaria in the world. Two of the most widely distributed vectors in Africa are Anopheles gambiae s.s. Giles and Anopheles arabiensis Patton, which are both members of a species complex of eight closely related, morphologically indistinguishable species known collectively as Anopheles gambiae s.l. [ 11]. In many regions the larval habitats of An. gambiae s.s. and An. arabiensis are similar, and in fact, the two species are often found within the same larval habitats [ 12- 14]. These larval habitats are generally smaller, temporary bodies of standing water persisting for about 20 to 40 days [ 12, 15], with rain being the main source of the water.
The locations of larval An. gambiae s.l. habitats are associated with certain environmental features of landscapes. Previous studies have found more larval habitats closer to streams [ 16] and in locations with agricultural land uses [ 16, 17]. Others have used a topographic wetness index (TWI) [ 18] to predict the locations of larval habitats, finding more larval habitats in locations having a combination of greater upslope area contributing to drainage and less slope [ 19- 21]. The influence of soil types on the presence of larval habitats has largely been ignored, although Bøgh and colleagues [ 22] found larval habitats exclusively in alluvial soils in The Gambia. Finally, seasonal differences in rainfall likely influence the number of larval habitats on the landscape [ 12, 16, 19, 20, 23].
The objectives of this study were to create a model for predicting larval An. gambiae s.l. habitat locations using landscape variables that predict the likelihood of standing water bodies, and to account for seasonal changes in habitat probability based on accumulated precipitation. A model for accurately predicting the locations of malaria vector larval habitats has multiple utilities. First, it allowed us to investigate the links between larval habitat distribution and adult malaria vector distribution across a large landscape where manually mapping the larval habitats is infeasible (McCann et al. in preparation). Additionally, such a model could be useful for malaria control programs, allowing program managers to focus their efforts to areas where larval habitats are most likely to occur.
The Asembo region of Rarieda District in western Kenya (Figure 1A) is a rural community of about 60,000 people covering about 200 km 2. Most of the residents are subsistence farmers, and the landscape is largely dominated by small-scale agriculture. Small plots of land generally surround family-based groups of houses, or compounds, further arranged into villages. While the compounds are highly dispersed within villages, the boundaries between villages are often discernable only by residents [ 24] (Figure 1B). Asembo sits in the lowlands along the shores of Lake Victoria, with elevations ranging from 1,100 m to 1,400 m above sea level and low topographic relief. Networks of streams run across the region and drain into Lake Victoria. Farmland is common in these low-lying drainage basins, as well as throughout the region. Houses are mostly absent within 100 m of the streams. Rainfall is seasonally bimodal but may occur year round, with monthly precipitation totals ranging from 7 to 490 mm and yearly totals ranging from 1,100 to 1,800 mm from 2003 through 2012.
Study site and survey locations. (A) Map of Kenya with red square indicating location of the study region in western Kenya. (B) Map showing the boundaries of Asembo and the streams within the community, with black dots representing all households. (C) 10 by 10 km study site shown as dashed black line with thirty-one 500 by 500 m quadrats surveyed for larval habitats shown as red boxes. (D) Location of the neighboring villages of Aduoyo-Miyare and Nguka, sites of the 15 monthly ground surveys for larval habitats, shown in red.
Malaria is holoendemic in Asembo, with parasitemia rates in children under 5 being around 50% in 2009 [ 25]. Similar to rainfall patterns, malaria transmission occurs year round, with seasonal peaks in May-July and October-November. The predominant species of malaria is P. falciparum. Two of the primary malaria vectors in the region are An. gambiae s.s. and An. arabiensis, the only two members of the An. gambiae s.l. species complex found here. The other primary malaria vector in the region is Anopheles funestus Giles. However, An. funestus and An. gambiae s.l. larvae do not generally occupy the same habitats, as the larval habitats of An. funestus are generally larger and more permanent than those of An. gambiae s.l. [ 12]. Larval An. gambiae s.l. habitats are numerous and widespread in Asembo yet heterogeneously distributed [ 16]. This makes it difficult to establish a relationship between larval habitats and the spatial distribution of the adult vectors or malaria prevalence in people. A 10 by 10 km study site was defined within Asembo to examine variation in the determinants of larval habitat location across a relatively large area. Because we wanted to include the lakeshore in the study site, the southern border fell largely within the lake, leaving 96.43 km 2 of actual landmass in the 10 by 10 km site.
The 10 by 10 km study site was divided into 500 by 500 m quadrats for larval An. gambiae s.l. habitat ground surveys. After excluding the quadrats that fell completely in the lake, we selected quadrats for larval habitat ground surveys using spatially stratified random sampling from the remaining 393 quadrats (Figure 1C). For spatial stratification the 10 by 10 km area was divided into 2 by 2 km blocks. The 500 by 500 m quadrats were randomly selected from groups defined by the 2 by 2 km blocks. Spatial stratification was implemented to avoid the problem of sampling a cluster of quadrats in a certain area of the grid, assuring spatial variation in the predictor variables [ 26]. The time required for surveying a quadrat varied greatly according to the number of larval habitats in the quadrat, which was not known a priori. Therefore, we surveyed as many quadrats as possible during the targeted time frame, which was the end of the long rainy season to coincide with the peak An. gambiae s.l. population level [ 27- 29]. Thus, 31 quadrats were surveyed exhaustively over 22 days between 17 May 2011 and 4 July 2011.
All potential larval An. gambiae s.l. habitats found in the quadrats were georeferenced with GPS units. Six field workers spaced 20 m from each other walked from one end of a quadrat to the other, using ArcPad (ESRI, Redlands, CA) on a GPS unit for navigating the borders of the quadrat. This was repeated until the entire quadrat was covered, usually in four to five passes. This approach allowed us to say, with certainty, where habitats were absent during the survey. In addition to recording the locations of each larval habitat, we recorded the presence or absence of Anopheles larvae. Larval An. gambiae s.l. habitats were defined as any standing body of water, regardless of whether Anopheles were present on the day of the ground survey, and falling under the following categories: drainage channel, burrow pit, rain pool, runoff, cluster of hoof prints, stream bed pool, pond/reservoir, wet meadow, well and tire track [ 15]. For a subset of habitats (the first five Anopheles-positive habitats for each of the four to five passes across each quadrat), Anopheles larvae and pupae were collected to confirm that the habitats were being used by An. gambiae s.l. All visible Anopheles larvae, up to a maximum of 20, were collected using a 300 ml dipper or plastic pipette as appropriate according to the size of the habitat. The specimens were transported to the lab for species identification. Larvae were raised to fourth-stage instars for identification, while pupae were allowed to eclose as adults before identification. All identifications were done according to Gillies and Coetzee [ 30].
To capture variation in habitat location across time due to seasonal rainfall patterns, additional ground surveys were conducted monthly in two neighboring villages, Aduoyo-Miyare and Nguka, covering 6.22 km 2 within the 10 by 10 km study site (Figure 1D). Two local field workers with extensive knowledge of the villages walked throughout the whole of each village over the course of one to three days, depending on the number of habitats encountered, each month from April 2011 through June 2012. Potential An. gambiae s.l. habitats were defined and recorded as above. Thus, we had the ability to say where habitats were present and absent within the two villages each month.
Spatial data for soils, land use-land cover (LULC), distance to the nearest stream, and TWI were created across the study site. These data were assembled in ArcGIS 10.0 (ESRI, Redlands, CA) in raster data structures with a spatial resolution of 20 m. All four datasets were treated as constant over time. Soil data were taken from the 1:1,000,000 exploratory soil map of Kenya, compiled by the Kenya Soil Survey in 1980 [ 31]. The three soil types in Asembo were 1) friable clay, 2) friable clay/sandy clay loam, and 3) firm, silty clay/clay. Of these soil types, friable clay drains more quickly, and firm, silty clay/clay drains more slowly. A satellite image from the IKONOS-2 sensor was used to create the LULC classification. Briefly, unsupervised classification was done using the K-means method [ 32] in ENVI 4.8 (Exelis Visual Information Solutions, Boulder, CO). Classes were combined into a binary data layer of agricultural or non-agricultural land use. All streams in Asembo were mapped using GPS units, and the Euclidean distance in meters to the nearest stream was calculated.
The TWI data were derived from a digital elevation model (DEM) of the study site. The DEM was created using local universal kriging to interpolate 11,130 GPS elevation records previously taken within Asembo [ 33, 34]. The ArcGIS extension TauDEM 5.0 (Tarboton, Utah State University) was used to calculate a TWI. First, slope and flow direction were calculated using the deterministic infinite-node algorithm recommended by Tarboton, which is robust yet easily implemented [ 35]. Contributing area catchments were calculated using the flow direction and slope data. Finally, TWI was calculated as the ratio of the slope to the contributing area, and the values were rescaled to the range by taking:
Because the TWI was calculated as the ratio of slope to contributing area, the lowest value (0) represented the wettest locations, while the highest TWI value (100) represented the driest areas.
Daily precipitation totals for March 2011 to July 2012, as measured by the weather station at the Kisumu Airport (about 40 km east of Asembo), were downloaded from the National Climatic Data Center’s Global Summary of Day (GSoD) database (Figures 2 and 3). For missing daily data at the Kisumu weather station (n = 19 of 489 days), the inverse distance weighted mean of surrounding GSoD weather stations (within 250 km) was used. Cumulative n-day precipitation totals were calculated by summing the precipitation total of a given day with the previous n days for n = 0 to 30 days. Each n-day precipitation total had a daily temporal resolution, treated as spatially constant across our study site. From the resulting 31 measures of cumulative precipitation, we selected the best cumulative n-day total for each model based on the criteria outlined below.
Precipitation data from the Kisumu airport for March 2011 through July 2012. Each point along the x-axis represents one day. The vertical bars show the daily precipitation total. The solid line shows the cumulative 30-day precipitation total (i.e. the level of the line for each day represents the summation of precipitation totals over the previous 30 days), which was used in the logistic regression model for the 15-month dataset. Similarly, the dashed line shows the daily fluctuation of the cumulative 21-day total used in the random forest model of the same data. Circles and squares show the last day of ground surveys for each month in Aduoyo-Miyare and Nguka. Kisumu airport is about 40 km east of Asembo.
Precipitation data from the Kisumu airport for 2 May 2011 to 4 July 2012. Each point along the x-axis represents one day. The vertical bars show the daily precipitation total. The solid line shows the cumulative 6-day precipitation total (i.e. the level of the line for each day represents the summation of precipitation totals over the previous 6 days), which was used in the logistic regression model of the 10 by 10 km dataset. Similarly, the dashed line shows the daily fluctuation of the cumulative 14-day precipitation total used in the random forest model of the same data. Squares and circles show the days of the ground surveys in the 31 quadrats. Kisumu airport is about 40 km east of Asembo.
We used two approaches for modeling the distribution of larval habitats across the landscape, logistic regression and random forest[ 36]. Logistic regression is commonly used in species distribution modeling [ 20, 21, 37]. Ecologists have recently started using the random forest method as well, because it does not require any assumptions about the distribution of the data [ 38, 39]. Random forest is a machine learning classification method that extends classification and regression tree (CART) approaches, which work by recursive binary partitioning of the data space into increasingly homogenous regions [ 39, 40]. Random forest works by fitting and combining many CARTs to create a more accurate prediction [ 36, 39].
Both methods were used separately on the two datasets (one from the 10 by 10 km area and the other from the 15 monthly ground surveys in Aduoyo-Miyare and Nguka). For both methods, the unit of analysis was a 20 m pixel. Because each of the 31 quadrats in the 10 by 10 km area was surveyed once, each pixel in that dataset had a single value for all of the variables described above (habitat presence/absence, TWI, soil, LULC, distance to stream, and 31 values of n-day precipitation). The value of n-day precipitation was based on the day a given quadrat was surveyed (Figure 3). In the Aduoyo-Miyare and Nguka dataset, each pixel was repeated 15 times. The value of n-day precipitation was based on the final day of ground surveys, which took only 1 to 3 days, each month (Figure 2). Habitat presence/absence was determined by the ground survey data for each month, while TWI, soil, LULC and distance to stream were constant for a given pixel across all 15 months.
We built a series of candidate logistic regression models to select the most useful predictor variables. To determine which n-day precipitation measure to use, each cumulative precipitation measure was used alone as the predictor variable in separate regression models. Of these, the model with the lowest BIC determined the cumulative precipitation measure used in the subsequent logistic regression candidate models. The five predictor variables (TWI, soil, LULC, distance to stream, and precipitation) were then used in all 31 possible combinations to build candidate logistic regression models for each of the two datasets. To restrict the candidate model sets to relatively simple models with easily interpretable parameters, interactions among the five predictor variables were not included. The top models were again selected according to the lowest BIC. While the locations of larval An. gambiae s.l. habitats are generally clustered, we did not account for spatial autocorrelation in the logistic regression models presented here. Previous studies modeling larval habitats have found similar results for logistic regression models with and without parameters accounting for spatial autocorrelation [ 19- 21]. The logistic regression models were implemented in the statistical software R 2.14.2 (R Development Core Team, Vienna, Austria).
We implemented the random forest approach using the R package ‘randomForest’ [ 41]. The best cumulative precipitation measure for use in the random forest models was determined by the mean decrease in the Gini impurity criterion when removing the variable from the full model (TWI, soil, LULC, distance to stream, and precipitation). The cumulative precipitation measure with the highest mean decrease in the Gini impurity criterion was used in the final random forest models. Gini impurity criteria were also used to measure variable importance in the random forest models. A greater mean decrease in the Gini impurity criterion suggests a stronger association with the response variable [ 36].
The top models from both approaches within each dataset were evaluated by determining their accuracy at predicting larval habitat presence and absence for holdout data. Fifty percent of each dataset was randomly selected as a holdout dataset before model building. Evaluation of model accuracy required the selection of a threshold at which to convert predicted probabilities into larval habitat presence or absence. Because threshold specific accuracy statistics can be sensitive to the threshold used for conversion, we generated an optimal threshold value by minimizing the absolute value of the difference between sensitivity and specificity [ 42]. This approach was chosen because both sensitivity and specificity were equally important for the intended application of the predictive models. To assess the performance of each model, we calculated the sensitivity, specificity, percent correctly classified (PCC), and kappa of each approach at each of the thresholds from the methods above. We also calculated the threshold-independent area under the curve (AUC) statistic [ 43] from receiver operating characteristics (ROC) plots using the R package ‘SDMTools’ [ 44].
Finally, we calculated Pearson’s correlation coefficient among the cumulative precipitation measures to assess differences among the temporal-resolution/modeling-approach combinations. To quantify the contribution of cumulative 30-day precipitation to variation in the number of habitats found each month in Aduoyo-Miyare and Nguka, we used simple linear regression.
In the 31 sampling quadrats selected from the 10 by 10 km study site, we recorded the locations of 1,673 larval An. gambiae s.l. habitats. Six of the quadrats did not have any larval habitats, while the mean number of habitats per 500 by 500 m quadrat was 54. Anopheles larvae were present in 921 of the 1,673 habitats on the day each habitat was recorded. As detailed in the methods, Anopheles larvae and pupae were collected from 141 of the habitats, 77% of which were occupied by An. gambiae s.l. on the day of collection. Most of the larvae and pupae (79%) were identified as An. gambiae s.l. The other species collected were An. funestus (1.1%), Anopheles coustani Laveran (6.7%), Anopheles rufipes (Gough) (5.3%), Anopheles maculipalpis Giles (2.5%) and Anopheles pharoensis/squamosus Theobald (3.9%).
In 15 monthly ground surveys in Aduoyo-Miyare and Nguka, a total of 6,770 larval An. gambiae s.l. habitats were recorded. The number of larval habitats in this area varied by month, ranging from 104 to 953 with a mean of 451. The number of larval habitats recorded in the two villages each month increased with increasing cumulative 30-day precipitation on the final day of ground surveys each month, though considerable variation was observed (R 2 = 0.1931, p = 0.1012; Figure 4).
Number of larval habitats by cumulative precipitation. Scatter plot of the number of habitats recorded in Aduoyo-Miyare and Nguka each month (n = 15) by the cumulative 30-day precipitation for the last day of ground surveys for that month. Line shows prediction of linear regression, R 2 = 0.1931, p = 0.1012. The red and blue boxes highlight variation in the residual error discussed further in the text.
The best cumulative precipitation total to use in the models differed between the datasets and between the modeling approaches. For the 15 monthly ground surveys in Aduoyo-Miyare and Nguka, the logistic regression model for 30-day cumulative precipitation had the lowest BIC within the precipitation candidate models (Figure 2), whereas the random forest model using the cumulative 21-day precipitation (Figure 2) had the highest mean decrease in the Gini impurity criterion. For the 10 by 10 km data, the logistic regression model for 6-day cumulative precipitation had the lowest BIC within the precipitation candidate models (Figure 3), while the random forest model using the cumulative 14-day precipitation (Figure 3) had the highest mean decrease in the Gini impurity criterion. Each of these precipitation measures was used within its respective dataset and modeling approach moving forward. However, it should be noted that these cumulative precipitation totals were moderately to highly correlated with each other (Table 1).
Pearson’s correlation coefficient (r) matrix of cumulative precipitation measures
| 6leftbottomA) 17 May - 4 July 2011 (temporal scale for 10 by 10 km dataset) | |||||
| leftbottom6-Day | leftbottom0.411 | leftbottom1.000 | leftbottom0.721 | leftbottom0.477 | leftbottom0.611 |
| leftbottom14-Day | leftbottom0.332 | leftbottom0.721 | leftbottom1.000 | leftbottom0.794 | leftbottom0.578 |
| 6leftbottomB) April 2011 - June 2012 (temporal scale for Aduoyo-Miyare and Nguka dataset) | |||||
| leftbottom21-Day | leftbottom0.408 | leftbottom0.555 | leftbottom0.841 | leftbottom1.000 | leftbottom0.979 |
| left30-Day | left0.336 | left0.537 | left0.819 | left0.979 | left1.000 |
(A) For the 22 days of larval habitat ground surveys in the 10 by 10 km area from 17 May to 4 July 2011, and (B) For the last day of larval habitat ground surveys each month in Aduoyo-Miyare and Nguka from April 2011 to June 2012. 0-day refers to precipitation total for the day of ground surveys, while 6-day, 14-day, 21-day, and 30-day refer to the cumulative precipitation total for the day of the ground surveys plus the previous 6 days, 14 days, 21 days, or 30 days, respectively.
The environmental variables used in the best logistic regression models for predicting the locations of larval An. gambiae s.l. habitats differed slightly between the datasets. For the 15-month dataset from Aduoyo-Miyare and Nguka, the logistic regression model with the lowest BIC included all five of the variables (Table 2). No other model had a ΔBIC less than 20. Larval habitats were more likely to be found in locations with a lower TWI (i.e. wetter because of a lower slope to contributing area ratio), closer to streams, in agricultural land use, and in the friable clay/sandy clay loam soil type (Table 3). The probability of a larval habitat increased with increasing cumulative 30-day precipitation (Table 3).
Top four logistic regression candidate models for the 15 monthly ground surveys in Aduoyo-Miyare and Nguka
| leftbottomTWI + LULC + DS + Soil + Precip. | leftbottom45933.1 | leftbottomNA | leftbottom0.9999 |
| leftbottomTWI + DS + Soil + Precip. | leftbottom45955.6 | leftbottom22.5 | leftbottom<0.001 |
| leftbottomTWI + LULC + DS + Precip. | leftbottom46045.5 | leftbottom112.3 | leftbottom<0.001 |
| leftTWI + DS + Precip. | left46073.7 | left140.5 | left<0.001 |
Based on Bayesian information criterion (BIC), in order of increasing difference in BIC from the top model (ΔBIC) and decreasing BIC weight, w. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; Precip., cumulative 30-day precipitation total; NA, not applicable.
Odds ratios for top logistic regression model from the Aduoyo-Miyare and Nguka data (15 monthly surveys)
| leftbottom(Intercept) | leftbottom0.0378 | leftbottom0.0334 | leftbottom0.0426 |
| leftbottomTWI | leftbottom0.9365 | leftbottom0.9324 | leftbottom0.9405 |
| leftbottomLULC, Ag:NonAg | leftbottom1.3371 | leftbottom1.2096 | leftbottom1.4780 |
| leftbottomDS | leftbottom0.9980 | leftbottom0.9978 | leftbottom0.9981 |
| leftbottomSoil, Type3:Type2 | leftbottom0.7127 | leftbottom0.6715 | leftbottom0.7564 |
| leftPrecip. | left1.0033 | left1.0029 | left1.0036 |
Odds ratios are presented with 95% CI. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; Precip., cumulative 30-day precipitation total; Ag:NonAg, odds ratio of agricultural to nonagricultural LULC; Type3:Type2, odds ratio of the firm, silty clay/clay soil type to the friable clay/sandy clay loam soil type.
For the 10 by 10 km dataset, the logistic regression model with the lowest BIC included four of the variables (TWI, distance to stream, soil type, and cumulative 6-day precipitation). No other model had a ΔBIC less than 9 (Table 4). Larval habitats were again more likely to be found in locations with a lower TWI, closer to streams, and in the friable clay/sandy clay loam soil type (Table 5). Counterintuitively, the probability of a larval habitat according to this model decreased with increasing cumulative 6-day precipitation.
The top five logistic regression candidate models for the 10 by 10 km area
| leftbottomTWI + DS + Soil + Precip. | leftbottom6482.6 | leftbottom0 | leftbottom0.9822 |
| leftbottomTWI + LULC + DS + Soil + Precip. | leftbottom6491.8 | leftbottom9.3 | leftbottom0.0096 |
| leftbottomTWI + DS + Precip. | leftbottom6492.2 | leftbottom9.6 | leftbottom0.0082 |
| leftbottomTWI + LULC + DS + Precip. | leftbottom6502.0 | leftbottom19.4 | leftbottom< 0.001 |
| leftTWI + DS + Soil | left6658.0 | left175.4 | left< 0.001 |
Based on Bayesian information criterion (BIC), in order of increasing difference in BIC from the top model (ΔBIC) and decreasing BIC weight, w. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; Precip., cumulative 30-day precipitation total.
Odds ratios for the top logistic regression model for the 10 by 10 km area data
| leftbottom(Intercept) | leftbottom0.3156 | leftbottom0.2472 | leftbottom0.4030 |
| leftbottomTWI | leftbottom0.9117 | leftbottom0.8988 | leftbottom0.9248 |
| leftbottomDS | leftbottom0.9973 | leftbottom0.9969 | leftbottom0.9977 |
| leftbottomSoil, Type2:Type1 | leftbottom1.9105 | leftbottom1.4994 | leftbottom2.4343 |
| leftbottomSoil, Type3:Type1 | leftbottom1.4970 | leftbottom1.2024 | leftbottom1.8639 |
| leftPrecip. | left0.9700 | left0.9656 | left0.9745 |
Odds ratios are presented with 95% CI. TWI, topographic wetness index; DS, distance to the nearest stream; Precip., cumulative 6-day precipitation total; Type2:Type1, odds ratio of the friable clay/sandy clay loam soil type to the friable clay soil type; Type3:Type1, odds ratio of the firm, silty clay/clay soil type to the friable clay soil type.
The most accurate model for predicting larval An. gambiae s.l. habitat locations in the 10 by 10 km area was the random forest method with all five variables (AUC = 0.864). In this model, TWI had the greatest mean decrease in the Gini impurity criterion, followed by distance to stream, precipitation, soil and LULC (Figure 5). Removing the cumulative 14-day precipitation from this model reduced the accuracy of the model (Table 6). The best logistic regression model for the 10 by 10 km area was less accurate than the random forest model when evaluated against the holdout data (Table 6). When the probabilities of larval habitat location across the entire 10 by 10 km study site were estimated using the most accurate models from each method, the random forest model clearly produced a more heterogeneous landscape at a fine scale than that produced by the logistic regression method (Figure 6).
Mean decrease in the Gini impurity criterion as a measure of variable importance. Shown for the random forest model in the 10 by 10 km dataset. Greater decreases indicate a stronger association with the response. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; Precip., cumulative 14-day precipitation total.
Comparison of models predicting the presence of larval habitats in the 10 by 10 km area
| leftbottomRF: TWI + DS + Soil + LULC + Precip. | leftbottom0.864 | leftbottom0.806 | leftbottom0.789 | leftbottom0.790 | leftbottom0.216 |
| leftbottomRF: TWI + DS + Soil + LULC | leftbottom0.808 | leftbottom0.750 | leftbottom0.725 | leftbottom0.726 | leftbottom0.145 |
| leftLR: TWI + DS + Soil + Precip. | left0.799 | left0.748 | left0.709 | left0.711 | left0.133 |
Two random forest (RF) models are shown, with and without Precip. (the cumulative 14-day precipitation total). The best logistic regression (LR) model is also shown. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; AUC, area under the receiver operating curve; PCC, percent correctly classified.
Probability of larval habitat presence across the 10 by 10 km study site. As predicted by (A) the logistic regression model using TWI, distance to stream and soil type; and (B) the random forest model using TWI, distance to stream, soil type and LULC.
The most accurate model for predicting larval An. gambiae s.l. habitat locations over the 15 monthly surveys in Aduoyo-Miyare and Nguka was the random forest method with all five variables (Table 7). As above, removing the cumulative 21-day precipitation from the model reduced its accuracy. The best logistic regression model for the 15 monthly surveys in Aduoyo-Miyare and Nguka was less accurate than the random forest model when evaluated against the holdout data (Table 7).
Comparison of models predicting the presence of larval habitats in Aduoyo-Miyare and Nguka
| leftbottomRF: TWI + DS + Soil + LULC + Precip. | leftbottom0.871 | leftbottom0.820 | leftbottom0.773 | leftbottom0.774 | leftbottom0.102 |
| leftbottomRF: TWI + DS + Soil + LULC | leftbottom0.827 | leftbottom0.659 | leftbottom0.936 | leftbottom0.930 | leftbottom0.268 |
| leftLR: TWI + DS + Soil + Precip. | left0.733 | left0.621 | left0.704 | left0.703 | left0.045 |
Two random forest (RF) models are shown with and without Precip. (the cumulative 14-day precipitation total). The best logistic regression (LR) model is also shown. TWI, topographic wetness index; LULC, land use-land cover; DS, distance to the nearest stream; AUC, area under the receiver operating curve; PCC, percent correctly classified.
The use of models to predict the distribution of species is common in ecology [ 37], and novel approaches to building these models such as random forest have become more widely available in recent years. We used two methods to predict the probability of larval An. gambiae s.l. habitat across the landscape and over time, and the random forest method produced more accurate models than the logistic regression method. This may be due to the differences in predicted heterogeneity of larval habitats at fine scales between the two methods, which can be compared across the entire 10 by 10 km study site (Figure 6). Predictions from the random forest model are more fragmented, showing a closer proximity of high-probability locations to low-probability locations, relative to the estimates from the logistic regression model. The general pattern is similar for the predictions of both models at broad scales (Figure 6). However, the fine scale heterogeneity in the random forest estimates more closely reflects the nature of actual larval habitat distribution on the ground, where larval An. gambiae s.l. habitats are distributed as many small patches rather than one continuous, large patch.
The most important landscape variables for predicting larval habitat presence in these models were TWI and distance to the nearest stream. In the 10 by 10 km random forest model, the mean decreases in the Gini impurity criteria of TWI and distance to the nearest stream were much larger than those of LULC and soil (Figure 5), indicating a stronger association with the prediction of habitat presence [ 36]. In general practice, high quality soil and LULC data can be difficult to acquire. Given limited resources, our data suggest it is possible to build reasonably accurate larval habitat models without these two landscape variables. Nonetheless, soil and LULC do show an association with habitat presence according to the logistic regression models presented here.
An important question in the application of predictive larval habitat models is whether models parameterized with data for habitat locations in one season are applicable to another season [ 20]. The creation of larval An. gambiae s.l. habitats (which are temporary, small bodies of standing water) depends on rainfall, which varies seasonally across the range of the species complex. One strategy to model differences between seasons is to account for variation in precipitation. In the random forest models, accumulated precipitation was less important than TWI and distance to the nearest stream, but it was a more important predictor variable than soil and LULC (Figure 5). Additionally, we found more larval habitats in months with more precipitation compared to the same area in months with less precipitation (Figure 4). Thus, including accumulated precipitation in our models improved the accuracy of larval habitat location predictions. These results should be interpreted with caution given the use of a single location about 40 km from the study site as the source of precipitation data. Daily precipitation totals can be spatially heterogeneous at that scale. Despite this limitation, it is clear from previous work that variation in precipitation influences larval An. gambiae s.l. habitats [ 12, 16, 19, 20]. However, the relationship may be more complex than it first appears. For example, it may not be linear. Rather, the number of larval habitats may increase monotonically with accumulated precipitation up to a threshold, after which more of the water on the landscape flows as surface sheet or channeled water, which is unsuitable aquatic habitat for An. gambiae s.l. larvae. Additionally, different habitat types may respond differently to increasing accumulated precipitation. Standing water forming in drainage channels and stream bed pools may be described better by a threshold relationship than the water filling burrow pits, hoof prints and tire tracks, because the former develops from channel and sheet water made stationary by diminished water flows, whereas the latter forms from water accumulating into various catchments not associated with channels. These additional factors may explain some of the uneven residual error seen in Figure 4, where the 4 months in the red box falling above the fitted regression line have more larval habitats with a lower accumulated precipitation relative to the 3 months in the blue box falling below the fitted regression line.
The n-day cumulative precipitation measure used for each modeling approach within each dataset was selected according to the criteria outlined in the methods to maximize the predictive power of each model. However, comparing across modeling approaches within each dataset, the cumulative precipitation measures were highly correlated (Table 1). Thus, the choice between 21-day cumulative precipitation and 30-day cumulative precipitation, for example, may be less important in general practice than using either measure instead of the daily precipitation total (referred to as 0-day in Table 1). Comparing across datasets, which differed in temporal scale, the 6-day cumulative precipitation and the 30-day cumulative precipitation are only moderately correlated. Their differences in terms of model fit (BIC) could reflect temporal differences in hydrology on this landscape, but it may also reflect a limitation of the 10 by 10 km data collection (see below).A counterintuitive result of this study was that the odds of larval habitat presence decreased with increasing cumulative 6-day precipitation using the best logistic regression model of the 10 by 10 km data. Most likely this reflects a limitation of the 10 by 10 km data collection rather than the true influence of precipitation on larval habitat presence, given the range of cumulative 6-day precipitation over the 49-day period (1.5 mm – 51.1 mm; Figure 3). The sampling strategy for those data was designed to capture variation in landscape variables over space. While precipitation varied among the days of the ground surveys, we were not able to capture that variation over the full range of values for the landscape variables. Instead, the effect of accumulated precipitation in this particular model may be an indication of some other property differing between the quadrats sampled on days of higher and lower accumulated precipitation. Alternatively, the temporal scale over which larval habitats respond to variation in accumulated precipitation may be closer to monthly than daily. That is, ground surveys conducted at monthly intervals in the same area may be more likely to be different than daily samples within a month in the same area. As noted above, this may also reflect the use of a single location as the source of all precipitation data.
In addition to the use of precipitation data from one location, there were other limitations to this study. First, we did not account for spatial autocorrelation in the logistic regression models. Doing so may have slightly increased the confidence intervals associated with the parameters of those models, but it is unlikely to have changed the model comparisons or accuracy evaluations presented here. Previous studies modeling An. gambiae s.l. larval habitat locations have found similar results for logistic regression models with and without parameters accounting for spatial autocorrelation [ 19- 21]. Second, there were additional variables we could have included in our analysis, such as a model-based wetness index (MWI) or normalized difference vegetation index (NDVI). MWI are similar to TWI, but MWI use simulations of distributed catchment models to account for differences between groundwater gradients and surface gradients, thereby creating more accurate topographic data [ 45]. We used TWI here because it has performed well in other models of Anopheles larval habitats [ 19- 21], and is easily implemented compared with MWI. While our models using TWI showed high accuracy, further studies comparing the use of MWI and TWI in larval habitat modeling are needed. NDVI has also been associated with the distribution of malaria [ 46, 47], although some studies have found contradicting results [ 17, 48, 49]. NDVI is an indirect measure of available moisture, but NDVI values are additionally influenced by vegetation type and phenology. Thus, we used accumulated precipitation as a measure of available moisture.
Finally, the models developed here exclusively used physical and environmental factors as predictor variables, but the formation of larval An. gambiae s.l. habitats also depends on human behavior. For example, landowners in Asembo create small drainage channels around fields. Standing water left behind in the channels creates habitats for An. gambiae s.l. larvae [ 15]. The locations of these drainage channels are often in low-lying agricultural areas, and therefore our models were able to predict the locations of most of the drainage channels. However, drainage channels are not found in 100% of low-lying agricultural areas, probably in part because of individual variation in landowner decision-making. Larval habitats formed from burrow pits and aggregations of hoof prints are also subject to variation in human behavior. While our models were able to correctly predict the locations of most of these habitats, interactions between the physical landscape and human behavior likely account for some of the locations identified incorrectly by the models.
The sampling designs of these two datasets allowed us to address two complementary goals. The monthly surveys in Aduoyo-Miyare and Nguka captured variation in precipitation across both dry and rainy seasons in the same landscape. This provided a stronger logical basis for inferences about the relationship between seasonal variation in precipitation and variation in the location and number of larval habitats. The small spatial extent of Aduoyo-Miyare and Nguka made monthly surveys more feasible, but it also limited the applicability of the model results across a larger area. Conversely, limiting the ground surveys of the 31 quadrats from the 10 by 10 km study site to one season likely impeded our ability to infer much about the effect of precipitation on these data. On the other hand, concentrating our sampling effort to increase replication across space in the 31 quadrats captured more variation in landscape variables, allowing us to apply the results of models based on these data to a larger area.
As a general application, the spatially stratified sampling strategy used in the 10 by 10 km site could serve as a framework for creating predictive larval habitat models for larval control. Targeted larval control is often cited as a useful application of predictive larval habitat models [ 20, 21], and we agree that there is potential for this application. For example, malaria control programs could identify areas suited to environmental management such as filling in burrow pits and engineering drainage channels to drain more completely. Additionally, allowing larvicide application crews to focus on areas with a higher probability of larval habitat presence would reduce the time, and therefore the cost, of larviciding. However, models fitted to data from a single geographic location may have limited generalizability [ 50]. Malaria control programs could overcome this limitation by using spatially stratified random samples, repeated across a variable landscape, to build models that are useful over larger areas.
The authors declare that they have no competing interests.
RSM, MNB, JMV, JEG and EDW designed the study and implemented the data collection. RSM, JPM and DWM analyzed the data and assessed the models. All authors participated in the preparation of the manuscript, and read and approved the final manuscript.
We thank George Olang’ and Maurice Ombok for logistical support; Richard Owerah, Jared Sudhe, Evans Owino, Peter Owera and Michael Nyonga for assistance with field work; Nicole Smith for assistance with DEM interpolation; Saul Daniel Ddumba and Nathan Moore for advice about the GSoD precipitation data; the staff from the KEMRI/CDC field station in Kisian for support in organizing field work and for assistance with mosquito identification; and the residents of Asembo for their cooperation during field surveys. We also thank four anonymous reviewers for constructive criticisms of an earlier version of the manuscript. This work is published with the permission of the Director of the Kenya Medical Research Institute. This study was supported by a National Science Foundation Ecology of Infectious Diseases grant (grant no. EF- 0723770) with additional support from the Rhodes Thompson Memorial Fellowship Fund.
The opinions expressed by the authors of this article do not necessarily reflect the opinions of the U.S. Centers for Disease Control and Prevention.


With no friction, it should be even easier to break the triangle than on ice.
Shoot the ball very carefully.
The triangle will break, and the particles will spread out evenly.
You just need to have patience.
(Psst! If you get bored, use this slider to speed up time.)


Experimenting with in-browser logistic regression
URL load url ↗ Demo Box id Load from Box (under dev)


Sjoerd Mullender
Robbert van Renesse
Without question, this C program is the most obfuscated C program that has ever been received! Like all great contest entries, they result in a change of rules for the following year. To prevent a flood of similar programs, we requested that programs be non machine specific.
This program was selected for the 1987 t-shirt collection.
NOTE: If your machine is not a Vax-11 or pdp-11, this program will not execute correctly. In later years, machine dependent code was discouraged.
The C startup routine (via crt0.o) transfers control to a location named main. In this case, main just happens to be in the data area. The array of shorts, which has been further obfuscated by use of different data types, just happens to form a meaningful set of PDP-11 and Vax instructions. The first word is a PDP-11 branch instruction that branches to the rest of the PDP code. On the Vax main is called with the calls instruction which uses the first word of the subroutine as a mask of registers to be saved. So on the Vax the first word can be anything. The real Vax code starts with the second word. This small program makes direct calls to the write() Unix system call to produce a message on the screen. Can you guess what is printed? We knew you couldn’t! :-)
Copyright © 1984, Landon Curt Noll. All Rights Reserved. Permission for personal, educational or non-profit use is granted provided this this copyright and notice are included in its entirety and remains unaltered. All other uses must receive prior permission in writing from both Landon Curt Noll and Larry Bassel.


Police should make first move to mend fences with public
by Ding WangOnce seen as Asia’s finest, the Hong Kong police force has seen its public image and credibility plummet to rock bottom amid months of anti-government protests.
There are complicated and inter-related socio-political factors that have led to the increasingly tense relations between the police and the public, and between the government and the public.
One is that the boundary between “two systems” under “one country” has become increasingly blurred.
Beijing’s continuing political interference in the city’s affairs has given rise to widespread resentment among the people over the erosion of freedoms and the rule of law, and this has fueled the ongoing resistance movement.
Meanwhile, the escalating violence of “valiant” protesters and police brutality have produced a chemical reaction that has resulted in the collapse of mutual trust between members of the public and law enforcement.
In order to rebuild that mutual trust, I believe that the party vested with public power (i.e., the police) should first behave strictly in accordance with the law.
As Philip Dykes, the chairman of the Hong Kong Bar Association, wrote in an open letter to his members on Oct. 16, “all government’s acts must comply with the law and be subject to effective review and legal scrutiny”.
“There is nothing more corrosive to the Rule of Law than violence inflicted on members of the public by apparently unaccountable state agents,” he stressed.
Moreover, in order to repair relations with the public, the party vested with public power should take the value of human life very seriously, and respect human dignity when enforcing the law.
Police officers must never prevent and hamper humanitarian and legal rescue efforts carried out by medics and lawyers during the crackdown on protesters and the subsequent arrests, or else public resentment at law enforcement would only run deeper and deeper in the days ahead.
This article appeared in the Hong Kong Economic Journal on Nov 14
Translation by Alan Lee
– Contact us at<EMAIL_ADDRESS>JC/CG


"Blood laser barrage"
- Tears are replaced with the ability to charge and fire a powerful laser that travels in a straight line across the room, dealing a lot of damage (damage * 13) to any enemies it comes into contact with.
- Does not directly affect your Damage stat but hits enemies 13 times per charge, spawning Guppy flies with every hit.
- Synergizes very well with a lot of items including Tammy's Head, Tiny Planet, Inner Eye (Triple Shot), Mutant Spider (Quad Shot) and many others.
- An item only found in the devil room pool (Note: The Fallen boss pulls items from the Devil Room Pool)


SOME MOTHER'S SON A film review by Cynthia A. Bussiere Copyright 1997 Cynthia A. Bussiere
Ireland Feature Film
[This review was first posted to CINEMA-L on January 16, 1997.]
Director: Terry George (directorial debut) Screenwriters: Terry George and Jim Sheridan Producers: Jim Sheridan, Arthur Lappin, Edward Burke Associate Producer: Helen Mirren Music: Bill Whelan
Principal Cast Kathleen Quigley: Helen Mirren Annie Higgins: Fionnula Flanagan Gerard Quigley: Aidan Gillen Frank Higgins: David O'Hara Bobby Sands: David Lynch Danny Boyle: Ciaran Hinds
In the late seventies, when (or shortly before) Margaret Thatcher took office in her first term as Prime Minister, the British government embarked upon an apparently harder line toward the IRA and its imprisoned volunteers than had previously been adhered to. In the words of the film, the government adopted a three-pronged policy of "isolation, criminalization, and demoralization." "Isolation" referred to a policy of geographical containment under which the British set up a system of road blocks and closures to hamper the movement of the IRA throughout Northern Ireland. "Criminalization" referred to a policy under which IRA prisoners were denied previously held "special category" status as political prisoners and, instead, were treated by the criminal justice system as criminal terrorists. "Demoralization" referred to a policy of breaking the will and the morale of imprisoned individual IRA volunteers and of the organization as a whole.
SOME MOTHER'S SON, a docudrama, focusses on the criminalization and demoralization prongs of that policy, although the story is set in motion by the isolation prong. As such, the film telescopes into the time frame of its story various events that actually occurred in Belfast in response to the British hard line. These events culminated, as the film culminates, in a prisoners' hunger strike in the Maze prison between March 1 and October 2, 1981.
Twenty-one prisoners went on the strike, and ten of them died. Of the ten, only one was sentenced to a life term. (Terry George states on the film's web site that, had they lived, all of the nine non-lifers would have been released by 1987.) Among the dead was Bobby Sands, beautifully played in the film by John Lynch. Sands was a prominent IRA member and a leader, within the prison, of the hunger strike. During his prison term, he was elected to Parliament. He was a member of Parliament when he died.
The strike didn't end until the mothers of the remaining striking prisoners exercised their powers as their sons' legal guardians and authorized medical intervention to save their sons' lives. By doing so, these women effectively took their sons off the hunger strike. Of course, the dilemma faced by each of these women was horrendous: By saving their sons' lives, they betrayed their sons' convictions and destroyed the cause for which their sons had fought so long and hard in the first place.
Thus, the film centers on the common dilemma of two mothers whose imprisoned sons have joined the hunger strike. On the one hand, there is Annie Higgins, played by Fionnula Flanagan. Annie is a farmer's wife and dedicated IRA supporter. By the time the film opens, one of her sons has already been killed by the British and her remaining son Frank, played by David O'Hara, is an IRA fugitive from justice. Her third child is her fabulous smart-mouthed, quick-tempered, red-headed daughter named Theresa who takes after Annie. On the other hand, there is Kathleen Quigley, played by Helen Mirren. Kathleen is a teacher in a Catholic girls' school. She is politically uninvolved to the point of naivete, but she hates the IRA for the violence for which it stands. Like Annie, she has three children, two sons and a daughter. Unbeknownst to her, her elder son Gerard, played by Aidan Gillen, is an IRA volunteer. Her younger son, as shown by a brief but illuminating moment in the film, may or may not grow up to emulate his older brother. Her older daughter, played by Geraldine O'Rawe, has her hands full just coping. Both mothers are the mainstays of their families: Annie's husband is quite elderly and physically unable to work. Kathleen is a widow.
The film follows the gradual leftward politicization of Kathleen after her son Gerard's arrest, trial, and conviction on a charge of attempted murder. Annie's son Frank, tried jointly with Gerard, is also convicted. Both young men are sent to the Maze to serve their sentences. The courtroom scene in which these events take place is highly sensational and quite moving in its brutal brevity.
When Gerard immediately joins the series of small ongoing protests conducted by all the IRA prisoners in his cell block, Mirren portrays Kathleen as having no choice but to resort to political action to pro- tect and ensure his welfare in prison once the authorities begin retaliating against the prisoners. In effect, his cause becomes her cause, and she takes it up lovingly and ungrudgingly. The film realistically charts her progressive politicization in a way that smoothly and efficiently advances the story. Along the way she engages in her own small rebellion against her Catholic school employer and its principal, the formidable and implacable Sister Bernardina (the nightmare nun of everybody's Catholic school education in the fifties and early sixties).
Concurrently with Kathleen's political progress, and equally real- istically, the film charts the progress of the budding friendship between Kathleen and Annie. The two women begin as openly hostile to one another. They then gradually warm up to each other, bond very closely, and then gradually come a bit apart again. This progression features, among other things, a hilarious car ride along a beach that is meant to be funny, despite the fact that it stands for more than it is. Structurally, these personal developments are woven into Kathleen's political development, because it is Annie, whose political convictions and activism are in place by the time the film begins, who involves Kathleen in specific political activities.
Eventually, the series of protests at the prison lead to the hunger strike, which both Kathleen's and Annie's sons join. From this moment forward, both women move into the center of a mael- strom of the competing interests of the three major players in these events: The IRA, the British government, and the Catholic Church. The IRA is represented by Danny Boyle, the head of Sinn Fein in Belfast. Well and strongly played by Ciaran Hinds, Boyle is a fictional composite character, not the least part of which is probably Gerry Adams. Boyle stubbornly and aggres- sively supports the prisoners' cause and serves as their repre- sentative in the outside world. The British government is represented by a ruthless counter-terrorism law enforcement official (whose scenes are often intercut with newsclips of Thatcher advocating her government's line), who is determined to break the backs of the prisoners and who, to that end, pursues a zero-tolerance policy and a stance of non-negotiation. Steering a middle ground between these extreme factions, yet supportive of the prisoners, is the Catholic Church. Initially, the Church is represented by the head of the Catholic Church in Ireland, the cardinal who comes to the Maze to say Mass for the prisoners. Horrified at the conditions he witnesses there, he immediately holds a press conference, to the extreme embarrassment and consternation of the government. However, more often, more prosaically, and less dramatically, the Church is represented by Father Daly, the local parish priest, who has the heart-breaking tasks of giving the last sacrament to dying hunger strikers and presiding at their funeral Masses.
Ultimately, each mother, Kathleen and Annie, holds the life of her son in her hands and must make the decision to let him live or watch him die.
All of the performances in the film are strong, accurate, and (with a few exceptions) beautifully understated. Mirren is fabulous. Her Kathleen has a touching diffidence, a surprising naivete, and great quiet courage. She also has flashes of humor that are a tribute not only to Mirren's acting, but also to George's direction and Sheridan's and George's writing. Flanagan, whose character undergoes no political changes throughout the film, has ostensibly less to work with. However, because the depth of her character's personal changes is arguably greater than that of Kathleen's, she alternates brilliantly between hardened political ideologue and blooming, softening private person. In a fabulous piece of acting, her physical appearance reflects these develop- ments by actually changing from within in front of the camera. Hinds, whose strong, quiet presence defined the character of Frederick Wentworth in PERSUASION, is an incredibly compelling screen presence here. His Danny Boyle also radiates strength, but of a steely, threatening, potentially explosive kind. As such, his character is as implacable as that of his British anti-terrorist counterpart -- and that of Sister Bernardina's, to boot. His char- acter is, unfortunately, rather underwritten, especially since it appears that at least one of the real-life characters of which his is a composite had a greater hand in the prison protests than the film conveys. Finally, the actors -- Aidan Gillen, John O'Hara, and John Lynch -- who portray the prisoners display a stunning talent that we rarely see in American actors of the same age. Furthermore, they all have strong film credits. Gillen was in CIRCLE OF FRIENDS, O'Hara was in BRAVEHEART, and Lynch was in IN THE NAME OF THE FATHER. In prison, with their uncut hair and unshaven beards, these young men have Christlike appearances. This imagery is further reinforced by background shots of anti-British posters depicting a young man crucified on the Union Jack.
Bill Whelan's folk-based Irish music is a strong and palpable presence in the film. Whelan, a record producer and theatre and television composer, is probably best known as the composer -- in its entirety -- of RIVERDANCE THE STAGE SHOW.
The film itself is extremely well thought out and therefore very well made. Like Spike Lee's DO THE RIGHT THING, the filmmakers and the actors know what they're talking about. Just as DO THE RIGHT THING conveyed a deadly accurate anatomy of a race riot, SOME MOTHER'S SON conveys a deadly accurate anatomy of a political hunger strike in which the lives of twenty-one young men turn on things like the interpretation of words such as "right" and "privilege" regarding what clothes to wear.
Almost nothing is wasted in this film. Therefore, the smallest occurrences and visuals, such as an apple swiped from a lunch or a member of Parliament smoking a cigarette, take on signifi- cance. The film opens and closes with a view of the Irish Sea. In between, and throughout the film, shots of feet on firm ground, such as marchers' feet in a candlelight vigil, are everywhere. The film also makes sparing but effective use of slow motion. Additionally, some of the editing, such as the intercutting of a portable missile attack with a close-up of girls' feet doing an Irish hard-shoe dance, is stunningly effective. Finally, the two principal crowd scenes -- the candlelight vigil for Bobby Sands and Sands' funeral -- are outstandingly well done. (The funeral alone comprised an entire second unit.) These scenes were not only exceptionally well-choreographed; they were also extraordinarily moving in a way that crowd scenes just usually are not.
Even though the film is essentially a work of fiction, all of the major events it depicts actually took place. These include the hunger strike itself, the various prison protests -- including the infamous "dirty protest" -- leading up to the strike, the cardinal's visit to the prison, the election of Sands to Pariliament, the vigils, the deaths, the funerals, and the secret deals, talks, and negotia- tions surrounding the hunger strike. Furthermore, while the characters in the film, with the exception of Sands and a member of Parliament, are fictitious, they are all composite characters based on real people. It appears that Belfast viewers may be able to recognize and identify the people on whom the characters are based.
The film was shot on location in Skerries, a small village north of Dublin. Skerries stood in for Belfast, where the principal events of the film take place. The Maze, together with its H-blocks where the prisoners were housed, was recreated in some Dublin warehouses.
The bottom line of this film is the interplay between the personal and the political in the lives of its characters. Sometimes the two are inseparable, sometimes they are in sync with each other, and sometimes they are at odds with each other. There is no question that this film has an explicit republican bias. However, there is also no question that the issues raised in the film are ultimately resolved on personal grounds.
The review above was posted to the
rec.arts.movies.reviews newsgroup (de.rec.film.kritiken for German reviews).
The Internet Movie Database accepts no responsibility for the contents of the
review and has no editorial control. Unless stated otherwise, the copyright
belongs to the author.
Please direct comments/criticisms of the review to relevant newsgroups.
Broken URLs inthe reviews are the responsibility of the author.
The formatting of the review is likely to differ from the original due
to ASCII to HTML conversion.
Related links: index of all rec.arts.movies.reviews reviews


One Interface To Rule Them All
Apache Libcloud is Python library that abstracts away differences among multiple cloud provider APIs.
Supports more than 30 providers such as
--- layout: default title: Apache Libcloud is a standard Python library that abstracts away differences among multiple cloud provider APIs description: Apache Libcloud is a Python library that abstracts away differences among multiple cloud provider APIs ---
Supports more than 30 providers such as
Latest stable version: 0.13.2
pip install apache-libcloud
Or download it from our servers and install it manually.
For more compute examples, see documentation.
For more DNS examples, see documentation.
{{ post.title }} {% endunless %} {% endfor %}
You can also subscribe and stay up to date using our RSS / Atom feed.
Users mailing list<EMAIL_ADDRESS>Developers mailing list<EMAIL_ADDRESS>IRC channel: #libcloud on Freenode


{{debug this}}
{{{class_description}}}
This class has no methods.
This class has no properties.
This class has no constructors.


| 017:053 Khan | : |
And say to My slaves (i.e. the true believers of Islamic Monotheism) that they should (only) say those words that are the best. (Because) Shaitan (Satan) verily, sows disagreements among them. Surely, Shaitan (Satan) is to man a plain enemy. |
| 017:053 Maulana | : |
And say to My servants that they speak what is best. Surely the devil sows dissensions among them. The devil is surely an open enemy to man. |
| 017:053 Pickthal | : |
Tell My bondmen to speak that which is kindlier. Lo! the devil soweth discord among them. Lo! the devil is for man an open foe. |
| 017:053 Rashad | : |
Tell My servants to treat each other in the best possible manner, for the devil will always try to drive a wedge among them. Surely, the devil is man's most ardent enemy. |
| 017:053 Sarwar | : |
(Muhammad), tell My servants to say what is best. Satan sows dissension among them; he is the sworn enemy of human beings. |
| 017:053 Shakir | : |
And say to My servants (that) they speak that which is best; surely the Shaitan sows dissensions among them; surely the Shaitan is an open enemy to man. |
| 017:053 Sherali | : |
And say to MY servants that they should always speak what is best. Surely, Satan, stirs up discord among them. Surely, Satan is an open enemy to man. |
| 017:053 Yusufali | : |
Say to My servants that they should (only) say those things that are best: for Satan doth sow dissensions among them: For Satan is to man an avowed enemy. |
Concepts/Themes Covered: أركان الإسلام > التوحيد > الكافرون > عداوة الكفار Pillars of Islam> Islamic> The Disbelievers> The Disbelievers' Hostility الإيمان > الغيب > الشيطان > سلوكه الشيطاني Faith> The Unseen World> Satan, the Devil> His Devilish Conduct الإيمان > الغيب > الشيطان > وسوسته Faith> The Unseen World> Satan, the Devil> Satanic Temptations الدعوة إلى الله > الحكمة في الدعوة > المجادلة بالتي هي أحسن The Call for Allah> The Wisdom in the Mission> Argument through the Optimum Ways العمل > العمل الصالح > البشاشة Action(Work)> The Good Deeds>Cheerfulness العمل > العمل الصالح > حسن السلوك Action(Work)> The Good Deeds>Good Conduct العمل > العمل الصالح > قول التي هي أحسن Action(Work)> The Good Deeds>To Talk in a Friendly Manner الإنسان والعلاقات الأخلاقية > الأخلاق الحميدة > البشاشة والوداعة Man and The Moral Relations> Good Morals>Cheerfulness and Humbleness الإنسان والعلاقات الأخلاقية > الأخلاق الحميدة > السلوك الحسن Man and The Moral Relations> Good Morals>Good Conduct الإنسان والعلاقات الأخلاقية > الأخلاق الحميدة > قول التي هي أحسن Man and The Moral Relations> Good Morals>Talking in a Friendly Manner
This tool uses eight English parallel translations for the Holy Quran imported from Mr. Clay Smith. You can also visit Qurany Concepts for browsing the Quran themes or topics. If you want to search the Holy Quran for keywords, please visit Qurany Keywords.


Tour of Rust
Table of Contents
Previous
The End
That's all for now. Stay tuned for new chapters. I hope you enjoy the journey ahead!