source('../../Data/get_species_data.R')
DF = get_species_data()
print(head(DF))

DF$logArea = log(DF$area)

m1 = glm(n_species ~ logArea, family=poisson, data=DF)
print(summary(m1))

m2 = glm(n_species ~ logArea + dist_to_nearest_neighbor, family=poisson, data=DF)
print(summary(m2))

m3 = glm(n_species ~ logArea + dist_to_santa_cruz, family=poisson, data=DF) ## seems "best"
print(summary(m3))

m4 = glm(n_species ~ logArea + n_species_on_adjacent, family=poisson, data=DF)
print(summary(m4))

m0 = glm(n_species ~ 1, family=poisson, data=DF)
step(m0, scope='n_species ~ logArea + dist_to_nearest_neighbor + dist_to_santa_cruz + n_species_on_adjacent', direction='forward', data=DF)

## Look at the residuals:
##
DF$mPrediction = predict(m3, type='response')
DF = DF[order(DF$n_species), ]
print(DF[, c('island', 'n_species', 'mPrediction')])