athowes / beyond-borders Goto Github PK

Evaluating the suitability of spatial adjacency for small-area estimation

Home Page: https://athowes.github.io/beyond-borders

License: MIT License

R 55.56% Stan 12.33% TeX 32.11%

random-effects-model small-area-estimation spatial-statistics

beyond-borders's Introduction

Hi! I'm Adam, a Bayesian statistician primarily working as a contractor with the Center for Forecasting and Outbreak Analytics at the CDC. Previously, I completed my PhD in statistics at Imperial College London. For more information, and my blog, see athowes.github.io. To give an idea about my work, see:

Repository	Description
`epidist`	An R package for estimating the time between events with interval censoring and right truncation
`thesis`	My thesis, "Bayesian spatio-temporal methods for small-area estimation of HIV indicators", supervised by Seth Flaxman (Oxford, CS) and Jeff Eaton (Imperial, DIDE; Harvard, CCDD)
`beyond-borders`	Evaluating the suitability of spatial adjacency for small-area estimation (`arealutils`)
`multi-agyw`	Estimating HIV risk group proportions with a multinomial spatio-temporal model (`multi.utils`, UNAIDS)
`naomi-aghq`	Improving approximate Bayesian inference methods for extended latent Gaussian models (`inf.utils`)
`exp-growth`	Time-series methods for disease-agnostic metagenomic environmental threat detection

I'm always looking to improve at writing clear, concise, efficient, reproducible code. Any suggestions, let me know! I predominantly program in R. Over time, I'd like to become multilinual by learning Python (and maybe Julia?): advice welcome.
I'm currently reading: Trees, maps, and theorems, The Effective Manager, and miscellaneous resources about technical management, see reading.
You may be interested in these resources I've found useful, or what I'm learning about via tutorial.

beyond-borders's People

Contributors

Stargazers

Watchers

Forkers

theorashid

beyond-borders's Issues

Use the newer assessments of coverage a la Säilynoja et al.

See discussion (and code) here: athowes/multi-agyw#23

Determine effect of geometry irregularity

#' For IK simulated data, what's the difference between Besag, FCK and FIK over the different geometries?
irregularity_effect <- df_rho %>%
  filter(sim_model == "IK", inf_model %in% c("Besag", "FCK", "FIK")) %>%
  group_by(geometry) %>%
  mutate(crps = 1000 * crps_mean,
         crps_min = min(crps),
         diff = crps - crps_min) %>%
  select(inf_model, diff)

# FCK versus best (FIK)
irregularity_effect %>% filter(inf_model == "FCK") %>% mutate_if(is.numeric, round, 1)

# Besag versus best (FIK)
irregularity_effect %>% filter(inf_model == "Besag") %>% mutate_if(is.numeric, round, 1)

Use `orderly` bundles to start running some of these reports on a cluster

Could try using a maths cluster, or otherwise something different to DIDEHPC since I've had trouble with that in the past
Though need to get cluster packages onto it (bsae) -- might be difficult for something other than DIDEHPC

Need to project before using hexagonal sampling in `cv_fik-inla`

[1] "Begin LOO cross-valdiation of CIV2017PHIA"
[ failed     ]  2022-02-20 11:00:05
[ elapsed    ]  Failed report in 0.2224343 secs
Error in st_poly_sample(x, size = size, ..., type = type, by_polygon = by_polygon) : 
  hexagonal sampling on geographic coordinates not supported; consider projecting first
In addition: Warning message:
 Error in st_poly_sample(x, size = size, ..., type = type, by_polygon = by_polygon) : 
  hexagonal sampling on geographic coordinates not supported; consider projecting first

Both operands should be units objects in `cv_fck-inla`

[1] "Begin LOO cross-valdiation of CIV2017PHIA"
[ failed     ]  2022-02-20 10:58:47
[ elapsed    ]  Failed report in 0.2870765 secs
Error in Ops.units(m, 1e-08) : 
  both operands of the expression should be "units" objects
In addition: Warning messages:
1: `progress_estimated()` was deprecated in dplyr 1.0.0. 
2: In st_centroid.sf(sf) :
  st_centroid assumes attributes are constant over geometries of x
3: In st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon) :
 
 Error in Ops.units(m, 1e-08) : 
  both operands of the expression should be "units" objects

Implement centroid and integrated kernel models in `R-INLA`

Is it done in Pointless spatial modelling, which has code?

Ideally all models should be implemented with the same inference algorithm to remove that from the comparison.

Move `areal-sim` repository to here

Move `sim_grid`, `sim_civ` and `sim_tex` to one report

Move commonly used functions to `bsae`

It's good practise to avoid duplicating code as much as possible

Shouldn't be generating posterior simulations so many times in model assessment

Simulate posterior, take one column, then discard, then simulate posterior, take another column and so on. Costly. Should find alternative way to do it (simulate posterior and use all the columns).

Add area file depends in make folder

Analyze lengthscale recovery

Need to run N_sim = 20 Stan models to get this starting. Running now.

Add table of inferential models to README

Move old `areal-kernel` repository over to here

Make the posterior assessment leaner

`approxfun` error when running `assess-rhos` for vignette geometries

[ elapsed    ]  Failed report in 1.758353 mins
Error in approxfun(marginal$df$x, marginal$df$y, yleft = 0, yright = 0) : 
  need at least two non-NA values to interpolate

Set the length-scale as oracle parameter

Show that `R-INLA` or `rstan` inbuilt CV measures are not accurate

Fit models with CV script and compare!

Square posterior marginal functions in `bsae` with posterior marginal functions here

Perform inference for vignette geometries

Add quantile and coverage assessment to HIV study

Move `areal-hiv` repository to here

Move any relevant things from the `scripts` folder of `bsae`

Create geometries (i) - (iv) as `sf`

For all but (ii) can use the polygon points
For (ii) use concentric circles e.g. https://gis.stackexchange.com/questions/322423/create-a-circle-in-r-and-plot-it-in-a-shapefile

Comment on messing up fixed effects? (Hodges & Reich, 2010)

https://www4.stat.ncsu.edu/~bjreich/papers/FixedEffectsYouLove.pdf

Start fitting the Stan models too

Might want to suppress rstan sampling output: it's getting in the way of my nice progress bar :(

Add Stan diagnostics to the `ck_stan` and `ik_stan` fitting reports

See how you've done it in elgm-inf bud.

Create annotated equations version of the Besag model

Useful for presentations, communication, etc.

Simulate from the vignette geometries

Assess both the vignette and realistic

Option to fix variance parameter

Jeff mentioned that it's unlikely that the variance parameter can be recovered in some of the vignette geometries (they have three areas). Add an option in bsae to fix this parameter to the correct value (oracle) then redo the inference with this option set here.

Some uncertainty here as to whether this benefits some of the methods more than others.

Add `vignette` and `realistic` toggle options for assessment scripts

Import both as depends using latest(vignette = TRUE) then use if statement in script.R

Investiate urban and rural for HIV

#' Largest five districts by population according to Wikipedia
big5 <- list(
  # https://en.wikipedia.org/wiki/Regions_of_Ivory_Coast
  "Cote d'Ivoire" = c("Abidjan", "Lagunes", "Woroba", "Savanes", "Lacs"),
  # https://en.wikipedia.org/wiki/Districts_of_Malawi
  "Malawi" = c("Lilongwe", "Blantyre", "Mzimba", "Mangochi", "Zomba"),
  # https://en.wikipedia.org/wiki/Districts_of_Tanzania
  "Tanzania" = c("Temeke", "Kinondoni", "Ilala", "Geita", "Sengerema"),
  # https://en.wikipedia.org/wiki/Districts_of_Zimbabwe
  "Zimbabwe" = c("Harare", "Shurugwi", "Bulawayo", "Mutare", "Gokwe South")
)

Look at the score cloropleths? Score by some measure of urbanicity?

Sort out the Tanzania islands issue from last time

Remove islands?

  nb <- sf_to_nb(tz) # Neighbourhood list
  comp <- spdep::n.comp.nb(nb) # Connected components
  table(comp$comp.id) # The number of districts in each component
  is_mainland <- comp$comp.id == 1
  tz_no_islands <- tz[is_mainland, ] # Only those on the mainland
  tz_no_islands$survey_id <- "TZ2012AIS-no-islands" # Change the ID so that fit_id works

Get DHS reports from naomi-orderly
They'll be more up to date, and shouldn't have to be making any alterations from the basic format