Malaria

import folium
import geopandas as gpd
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.colors import rgb2hex

from pyinla.model import *
from pyinla.raster import *
from pyinla.spde import *
from pyinla.utils import *

gambia = ro.r('library("geoR"); data(Gambia); gambia')

R[write to console]: --------------------------------------------------------------
 Analysis of Geostatistical Data
 For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
 geoR version 1.9-2 (built on 2022-08-09) is now loaded
--------------------------------------------------------------

agg = gambia.groupby(["x", "y"])

total = agg["pos"].size()
positive = agg["pos"].sum()
prevalence = positive / total

x, y = list(zip(*agg.groups.keys()))

d = dict(
    x=np.array(x),
    y=np.array(y),
    total=total.values,
    positive=positive.values,
    prevalence=prevalence.values,
)

gdf = gpd.GeoDataFrame.from_dict(
    d, geometry=gpd.points_from_xy(x=x, y=y), crs="+proj=utm +zone=28"
).to_crs("+proj=longlat +datum=WGS84")
geom = gdf.geometry.to_crs("+proj=longlat +datum=WGS84")
gdf["long"] = geom.centroid.x
gdf["lat"] = geom.centroid.y

/var/folders/l1/xkqt8tg56t325lgp9ctmslwm0000gn/T/ipykernel_92607/3175676654.py:5: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  gdf["long"] = geom.centroid.x
/var/folders/l1/xkqt8tg56t325lgp9ctmslwm0000gn/T/ipykernel_92607/3175676654.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  gdf["lat"] = geom.centroid.y

m = folium.Map(location=(gdf.lat[0], gdf.long[0]))
for i in range(len(gdf)):
    folium.Circle(radius=50, location=(gdf.lat[i], gdf.long[i])).add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook

raster = ro.r(
    'library(raster); getData(name = "alt", country = "GMB", mask = TRUE, path="./data")'
)
alt = raster_extract(raster, gdf[["long", "lat"]])
gdf["alt"] = alt

coo = gdf[["long", "lat"]]
mesh = mesh_2d(coo, max_edge=[0.1, 5], cutoff=0.01)

mesh.plot(figsize=(7, 5))

<AxesSubplot: >

ra = raster_aggregate(raster, fact=5, fun="mean")
dp = raster_to_points(ra)
coop = dp[:, :2]
Ap = make_projection_matrix(mesh=mesh, loc=coop)

A = make_projection_matrix(mesh=mesh, loc=coo)

spde = spde2_matern(mesh)
indexs = spde.make_index("s")

stk_e = inla_stack(
    tag="est",
    data=dict(y=gdf.positive, numtrials=gdf.total),
    A=A,
    effects=dict(b0=1, altitude=gdf.alt),
    s=indexs,
)
stk_p = inla_stack(
    tag="pred",
    data=dict(y=np.nan, numtrials=np.nan),
    A=Ap,
    effects=dict(b0=1, altitude=dp[:, 2]),
    s=indexs,
)
stk_full = combine_stacks(stk_e, stk_p)

formula = "y ~ 0 + b0 + altitude + f(s, model = spde)"

res = inla(
    formula,
    data=convert_r2py(stack_data(stk_full)) | dict(spde=spde.spde),
    family="binomial",
    n_trials="numtrials",
    control_family=dict(link="logit"),
    control_predictor=dict(compute=True, link=1, A=stack_A(stk_full)),
).improve_hyperpar()

res

Time used:
     = 1.78,  = 0.49,  = 0.0286,  = 2.3 
Fixed effects:
           mean   sd 0.025quant 0.5quant 0.975quant   mode kld
b0       -0.190 0.49     -1.048   -0.227      0.911 -0.278   0
altitude -0.012 0.01     -0.032   -0.012      0.009 -0.012   0

Random effects:
  Name	  Model
    s SPDE2 model

Model hyperparameters:
              mean    sd 0.025quant 0.5quant 0.975quant  mode
Theta1 for s -4.00 0.286      -4.56    -4.00      -3.43 -4.00
Theta2 for s  2.63 0.365       1.91     2.64       3.35  2.65

Marginal log-Likelihood:  -216.11 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

pred_index = stack_index(stk_full, tag="pred")

prevalence_summary = pd.DataFrame(res.get_summary("fitted.values"))
prevalence_summary = prevalence_summary.iloc[pred_index]

m = folium.Map(location=coop[0][::-1])
for i in range(len(coop)):
    cc = rgb2hex(cm.viridis(prevalence_summary["mean"].iloc[i]))
    folium.Circle(
        radius=2000, location=tuple(coop[i][::-1]), color=cc, fill=True, fill_color=cc
    ).add_to(m)
m

Make this Notebook Trusted to load map: File -> Trust Notebook