source('read_toxoplasmosis_data.R')

toxo = read_toxoplasmosis_data()

resp = cbind(toxo$Cases, toxo$Number_tested - toxo$Cases)
Rainfall = toxo$Rainfall / 1000

## Fit a linear model of rainfall:
##
m1 = glm(resp ~ Rainfall, family=binomial)
print(summary(m1))

## Look if number tested and Rainfall are correlated:
##
ct = cor.test(toxo$Number_tested, Rainfall)
print(ct)

## Fit a polynomial model of rainfall:
##
mp = glm(resp ~ Rainfall + I(Rainfall^2) + I(Rainfall^3), family=binomial)
print(summary(mp))

## Plot the models above:
##
plot(Rainfall, resp[, 1] / (resp[, 1] + resp[, 2]), xlab='rainfall/1000', ylab='probability of toxoplasmosis')
x_grid = seq(min(Rainfall), max(Rainfall), length.out=100)
m1_hat = predict(m1, newdata=data.frame(Rainfall=x_grid), type='response')
lines(x_grid, m1_hat, col='red')
mp_hat = predict(mp, newdata=data.frame(Rainfall=x_grid), type='response')
lines(x_grid, mp_hat, col='black')
grid()