CHWplacement
A package to compute optimal CHW geographical placement scenarios
CHWplacement-vignette.RmdThis R package contains functions that can be used to compute CHW placement scenarios, based on gridded population surfaces, friction surfaces and additional datasets such as GPS coordinates of existing health facilities.
These functions are based on the CLSCP optimisation programme by Current and Storbeck (1988). Further details on the methodology can be found in the associated manuscript (Champagne et al. (2022)).
The CHW placement scenarios are designed to find the minimum number of CHWs required to reach full population coverage under two constraints:
the travel time between the CHWs and the households remains below a given time
the population assigned to each CHW remains below a given threshold
The algorithm provides the total number of CHWs required and their GPS positions.
How to use the package
The package can be installed using the devtools package:
library(devtools)
install_github("SwissTPH/CHWplacement") First steps and preparation of model inputs
First the package needs to be loaded, with some additional libraries for visualisation.
library(CHWplacement)
library(rgeos)
#> Loading required package: sp
#> rgeos version: 0.5-9, (SVN revision 684)
#> GEOS runtime version: 3.9.1-CAPI-1.14.2
#> Please note that rgeos will be retired by the end of 2023,
#> plan transition to sf functions using GEOS at your earliest convenience.
#> GEOS using OverlayNG
#> Linking to sp version: 1.4-7
#> Polygon checking: TRUE
library(rasterVis)
#> Loading required package: lattice
library(ggplot2)
library(cowplot)The optimisation requires 5 inputs:
population.raster: a population gridded surface, representing population distribution
friction.raster: a friction surface (for example, the surface by Weiss et al. (2020) that can be downloaded from the malariaAtlas R package (Pfeffer et al. (2018))), representing the difficulty to cross one pixel in the maps (e.g. in minutes to cross one meter)
access.raster: a raster surface indicating areas to be included in the algorithm, e.g. accessibility to health facilities (for example, the surface by Weiss et al. (2020) that can be downloaded from the malariaAtlas R package (Pfeffer et al. (2018)))
popurb.raster: a raster surface indicating the positioning of urban areas (for example, as derived using the DefineUrban function)
shp: a polygon shapefile indicating the area to be considered for optimisation
For this example, we will demonstrate with fictitious data, but the user’s data should be used instead.
newproj= "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"
pop.dummy=raster::raster( matrix(c(rep(30,125), rep(90,100)), ncol=15))
fric.dummy=raster::raster(matrix(rep(0.02,225), ncol=15))
acc.dummy=raster::raster(matrix(c(rep(30,175), rep(90,50)), byrow=T, ncol=15))
shp.dummy <- as(raster::extent(pop.dummy), 'SpatialPolygons')
raster::crs(pop.dummy)=newproj
raster::crs(fric.dummy)=newproj
raster::crs(acc.dummy)=newproj
names(acc.dummy)=""
The raster file indicating the position of urban areas can be created following the “degree of urbanization approach” described in the 2017 World Bank report (The World Bank (2017)) to define urban and rural areas, with the DefineUrban function. The parameter rururb_cutoff indicates the minimum population density for each pixel considered as urban. The parameter min_urbsize indicates the minimum population in the total area of contiguous selected pixels to be considered as urban. The user can also use any raster file created with another methodology.
urb.dummy = DefineUrban( pop.dummy, rururb_cutoff = 50, min_urbsize = 1000)
#> Loading required namespace: igraph
#> [1] "Proportion of the population in urban areas is 70.59%"We can first visualise this fictitious data, using the following plotting packages:
plot_pop = gplot(pop.dummy)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +coord_equal()+scale_fill_manual(values=c("#FF6600", "#990000"), name = "Population")
plot_acc = gplot(acc.dummy)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +coord_equal()+scale_fill_manual(values=c("forestgreen", "darkorange"), name="Access")
plot_urb = gplot(urb.dummy)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +coord_equal()+scale_fill_manual( values=c("#990000"), name="Urban\npopulation")
plot_fric = gplot(fric.dummy)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +coord_equal()+scale_fill_manual(values=c("#FFCC00"), name="Friction")In this fictitious population, we have 12750 people, 9000 of them living in an urban area. The southern areas is considered to be hard-to-reach (orange, as opposed to green). This represents 720 individuals in the rural area, and 2340 in the urban area. The friction is uniform.
plot_grid(plot_pop, plot_fric, plot_acc, plot_urb)
Creating a CHW placement scenario on a simple example
These elements are combined to provide the minimal number of CHWs and their optimal geographical positioning to guarantee full population coverage verifying the following constraints:
Maximum travel time: the travel time between the CHW and their assigned inhabitants remains below a certain time (called “radius”)
Maximum population: the maximal number of inhabitant per CHWs remains below a certain capacity threshold, that can differ in urban (“max.treat.per.CHW.urban”) and rural (“max.treat.per.CHW.rural”) areas
Area selection: the analysis can be restricted to certain areas in the access.raster map via the “buffer” and “is.inside” parameters: for example if access.raster represents the travel time to health facilities in min, the analysis can restricted to areas situated at more than 60 min from any health facility (buffer=60 and is.inside=FALSE), or areas situated at less than 60 min from any health facility (buffer=60 and is.inside=TRUE)
The CHW plan can be defined with the following function:
CreateCHWplacement( population.raster=pop.dummy, friction.raster=fric.dummy,shp=shp.dummy,
name="temp", buffer=60, radius=600, capacity.name="",
access.raster=acc.dummy,popurb.raster=urb.dummy,
max.treat.per.CHW.urban=2000, max.treat.per.CHW.rural=1000,
max.CHW.per.pixel=10, is.inside = FALSE, filepath=".")
#> Warning in raster::projectRaster(friction.map, pop.map): input and ouput crs are
#> the same
#> Warning in raster::projectRaster(access.map, pop.map): input and ouput crs are
#> the same
#> as(<dsCMatrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(as(., "generalMatrix"), "TsparseMatrix") instead
#> Model name:
#> a linear program with 2550 decision variables and 100 constraintsThis function will create a file called temp.mps, stored in the repository indicated by filepath. This files contains the description of the optimisation problem, that can then be solved by an optimisation algorithm. For example, the problem can be solved in R with the lpSolveAPI package, with the following function:
CHWanswer=solve_CHWplacement(population.raster=pop.dummy, friction.raster=fric.dummy,shp=shp.dummy,
name="temp", buffer=60, radius=600, capacity.name="",
access.raster=acc.dummy,popurb.raster=urb.dummy,
is.inside = FALSE, filepath=".")
#> Warning in raster::projectRaster(friction.map, pop.map): input and ouput crs are
#> the same
#> Warning in raster::projectRaster(access.map, pop.map): input and ouput crs are
#> the same
#> [1] 33 CHWs are needed to cover the area. The package returns a new data.frame with one line per assigned CHW, with the associated GPS coordinates (x and y), whether their are positioned in rural areas (is.rural=1 in rural areas and 0 in urban areas), and their associated number of people assigned (capacity). The sum of all CHW capacities is above the total population in the area, due to overlapping catchment areas in this example.
CHWanswer
#> x y is.rural capacity
#> 1 0.7000000 0.23333333 0 1770
#> 2 0.8333333 0.23333333 0 900
#> 3 0.2333333 0.03333333 1 660We can represent these results graphically. To improve interpretability, we will plot a masked version of the population raster, where all areas not included in the analysis (following the buffer parameter) are replaced by NA. Such a map can be obtained as one of the outcomes of the function PrepareRasterFiles.
all.rasters=PrepareRasterFiles(population.raster=pop.dummy, friction.raster=fric.dummy,shp=shp.dummy,
buffer=60,
access.raster=acc.dummy,popurb.raster=urb.dummy,
is.inside = FALSE)
#> Warning in raster::projectRaster(friction.map, pop.map): input and ouput crs are
#> the same
#> Warning in raster::projectRaster(access.map, pop.map): input and ouput crs are
#> the sameWe use the rasters to represent the results graphically:
gplot(all.rasters$pop.map)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +scale_fill_manual(values=c("#FF6600", "#990000"), name="Population") +
geom_point(data=CHWanswer, aes(x=x, y=y))+
coord_equal()
We can also plot the catchment area of each CHW and see that all
population points are covered. These catchment areas correspond to a
distance of 600 min (radius) from the CHW location. The catchment areas
do overlap, and in practice, they would need to be adapted individually
knowing the local geographical context when CHWs are recruited.
footprint1=GetFootprintWalkingTime(pos=CHWanswer[1,1:2], radius=600, pop.map=all.rasters$pop.map, T.GC=all.rasters$T.GC)
footprint1b=raster::xyFromCell(all.rasters$pop.map,footprint1)
footprint2=GetFootprintWalkingTime(pos=CHWanswer[2,1:2], radius=600, pop.map=all.rasters$pop.map, T.GC=all.rasters$T.GC)
footprint2b=raster::xyFromCell(all.rasters$pop.map,footprint2)
footprint3=GetFootprintWalkingTime(pos=CHWanswer[3,1:2], radius=600, pop.map=all.rasters$pop.map, T.GC=all.rasters$T.GC)
footprint3b=raster::xyFromCell(all.rasters$pop.map,footprint3)
gplot(all.rasters$pop.map)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +scale_fill_manual(values=c("#FF6600", "#990000"), name="Population") +
geom_point(data=as.data.frame(footprint1b), aes(x=x, y=y), color="darkblue")+
geom_point(data=as.data.frame(footprint2b), aes(x=x, y=y), color="dodgerblue")+
geom_point(data=as.data.frame(footprint3b), aes(x=x, y=y), color="lightblue")+
geom_point(data=CHWanswer, aes(x=x, y=y), size=2)+
coord_equal()
Using the package in practice with large maps
Solving the model with the R software and the lpSolveAPI package can become very computationally intense for the problems to be solved in practice (larger maps and more difficult optimisation constraints). In this context, other software can be used, such as Gurobi Optimization, LLC (2020) or SCIP (Gamrath et al. (2020a), Gamrath et al. (2020b)).
Some limitations
Time limitations and size limitation are caused by the creation and .mps by the lpSolveAPI package. For large problems the .mps file creation can be lengthy. Additonally, there can’t be more than 9 999 999 decision variables, because the name of columns i.e. decision variables in .mps have a limited length (after that threshold, C10000001 and C10000002 become C1000000 which raises error).
Using Gurobi
Gurobi is a commercial whose performance is greater than lpSolveAPI on large problems. The CHW plan can be defined similarly to the previous examples, with the same function:
CreateCHWplacement( population.raster=pop.dummy,
friction.raster=fric.dummy,shp=shp.dummy,
name="Gurobitemp", buffer=60, radius=600, capacity.name="",
access.raster=acc.dummy,popurb.raster=urb.dummy,
max.treat.per.CHW.urban=2000, max.treat.per.CHW.rural=1000,
max.CHW.per.pixel=10, is.inside = FALSE, filepath=".")This function will create a file .mps, stored in the repository indicated by filepath. This files contains the description of the optimization problem for Gurobi.
It can be solved by Gurobi following Gurobi’s documentation. For example, in the command line, the user can use the following command:
gurobi_cl MIPGapAbs=3 ResultFile=mypath/temp_solution.sol mypath/myfile.mpsIn this example, the option MIPGapAbs enables the user to change the error tolerance for the algorithm convergence (“Gurobi gap”) and hence balance the precision of the solution with computational time. More details and other options can be found in Gurobi’s documentation.
After running Gurobi, a file temp_solution.sol should be created in the chosen folder and the next step is to decode it. Make sure you use exactly the same input rasters as the ones used in CreateCHWplacement.
CHWanswer = read_CHWplacement_sol( population.raster = pop.dummy,
friction.raster = fric.dummy,
access.raster = acc.dummy,
popurb.raster = urb.dummy,
shp = shp.dummy,
name = "Gurobitemp" ,
filepath = ".",
buffer=60, is.inside = F, radius=600, capacity.name = "",
write = F # creates a CSV file of the results
)This will return a dataframe file containing the CHW positions (it can be saved as .csv file with the option write=T).
Using SCIP
SCIP is an open source solver whose performance is greater than lpSolveAPI on large problems. It needs to be downloaded and installed here. You can indicate which version of SCIP was downloaded as follows:
SCIP_version = "7.0.3"The CHW plan can be defined similarly to the previous examples, with the following function:
# mps_directory_H2 = file.path("clscp_600_buffer60_capa")
CreateCHWplacement_SCIP( population.raster = pop.dummy,
friction.raster = fric.dummy,
access.raster = acc.dummy,
popurb.raster = urb.dummy,
shp = shp.dummy,
filepath = ".",
name = "SCIPtemp",
buffer = 60, radius = 600, capacity.name="",
max.treat.per.CHW.urban = 2000, max.treat.per.CHW.rural = 1000,
max.CHW.per.pixel = 10, is.inside = FALSE,
dualitygap = 0,
relativegap = F,
SCIP_version = SCIP_version)
#> Warning in raster::projectRaster(friction.map, pop.map): input and ouput crs are
#> the same
#> Warning in raster::projectRaster(access.map, pop.map): input and ouput crs are
#> the same
#> Model name:
#> a linear program with 2550 decision variables and 100 constraintsThis function will create a file called SCIPtemp.mps, stored in the repository indicated by filepath. This files contains the description of the optimization problem.
Then it will create a .txt file with the basic script you need in order to calculate the solution by using SCIP shell. The two additional options (dualitygap and relativegap are specific SCIP parameters as detailed in the function documentation). These parameters enable the user to change the error tolerance for the algorithm convergence (“Gurobi gap”) and hence balance the precision of the solution with computational time.
cat(readLines('clscp_600_buffer60_capa/SCIP_script_SCIPtemp.txt'), sep = '\n')
***** Manual use of SCIP version 7.0.3*****
** Please enter the lines below one by one on the SCIP shell **
set limits absgap 0
read ./clscp_600_buffer60_capa/SCIPtemp.mps
optimize
write solution ./clscp_600_buffer60_capa/SCIPtemp_solutions.sol
quit
***** Automated use of SCIP version 7.0.3( for Windows users ) *****
** Please enter the script below in windows powershell **
Start-Process -FilePath 'C:/Program Files/SCIPOptSuite 7.0.3/bin/scip.exe'
$wshell = New-Object -ComObject wscript.shell
$wshell.SendKeys('%{TAB}')
Sleep 2
$wshell.SendKeys( 'set limits absgap 0' )
$wshell.SendKeys('~')
Sleep .3
$wshell.SendKeys( 'read ./clscp_600_buffer60_capa/SCIPtemp.mps' )
$wshell.SendKeys('~')
Sleep 1
$wshell.SendKeys('optimize' )
$wshell.SendKeys('~')
Wait-Process -Id $wshell.id
$wshell.SendKeys( 'write solution ./clscp_600_buffer60_capa/SCIPtemp_solutions.sol' )
$wshell.SendKeys('~')
Sleep .5
$wshell.SendKeys('quit')
$wshell.SendKeys('~')This text file contains two parts. The first part (Manual use of SCIP) gives the commands to type in the SCIP command line to solve the problem. The second part (Automated use of SCIP version 7.0.3( for Windows users )) gives a script that can be directly copy-pasted in the Windows powershell to solve the optimisation problem with SCIP. This will create a .sol file in the same repository, containing the solution of the optimisation problem.
If you already have a .mps file, you can use the function CreateScriptSCIP, which only creates the .txt file. It can be useful if you want to recreate scripts for an older .mps file.
CreateScriptSCIP( mps_directory = "clscp_600_buffer60_capa", name = "SCIPtemp", dualitygap = 0, relativegap = F )Now you have a file temp_solution.sol and the next step is to decode it. Make sure you use exactly the same input rasters as the ones used in CreateCHWplacement_SCIP.
CHWanswer = read_CHWplacement_sol_SCIP( population.raster = pop.dummy,
friction.raster = fric.dummy,
access.raster = acc.dummy,
popurb.raster = urb.dummy,
shp = shp.dummy,
name = "SCIPtemp" ,
filepath = ".",
buffer=60, is.inside = F, radius=600, capacity.name = "",
write = F # creates a CSV file of the results
)
#> Warning in raster::projectRaster(friction.map, pop.map): input and ouput crs are
#> the same
#> Warning in raster::projectRaster(access.map, pop.map): input and ouput crs are
#> the sameThis will return a dataframe file containing the CHW positions (it can be saved as .csv file with the option write=T).
3 CHWs are needed to cover the area: the solution is slightly different from the previous one, but it equally verifies the chosen constraints (walking time and capacities).
CHWanswer
#> x y is.rural capacity
#> 1 0.7000000 0.03333333 0 1290
#> 2 0.2333333 0.03333333 1 600
#> 3 0.7666667 0.10000000 0 1530
gplot(all.rasters$pop.map)+
geom_tile(aes(fill=factor(value)), alpha=0.6) +scale_fill_manual(values=c("#FF6600", "#990000"), name="Population") +
geom_point(data=CHWanswer, aes(x=x, y=y))+
coord_equal()