Trends and drivers of cardiovascular health in coastal communities

[mnky9800n]

Here starts the discussion for the project "Trends and drivers of cardiovascular health in coastal communities"

[sebsfox]

MSOA21 admissions to hospital for cardiovascular disease.
MSOA21 shape file

Code to read sf file in for MSOA21

msoa21url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2021_Boundaries_Generalised_Clipped_EW_BGC_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"
msoa21_shp <- sf::read_sf(msoa21url)

[sebsfox]

Data dictionary

Data dictionary (vary by year and area):

• AreaCode – MSOA21CD code
• FYEAR – financial year
• CVD – number of hospital admissions where the patient’s most recent admission was related to a cardiovascular condition (ICD10 code)
• Resp – ditto, respiratory
• MSK – ditto, musculoskeletal
• Psy – ditto, psychological/behavioural disorders
• Population – population of all ages, based on Census 2011 populations
• Female - % female
• Male - % male
• Pop_working - % 18 to 65
• Pop18 - % 18 and below
• Pop65 - % 65 and above
• IMD_1 - % areas that are in the least deprived decile (might want to change this variable as many nulls)

Data dictionary (vary by area, mostly Census 2021 variables)

• Life_exp_m – male life expectancy
• Life_exp_f – female life expectancy
• No_qual - % with no qualifications
• Higher_qual - % with degree or higher
• Econ_inactive - % economically inactive
• Bad_health - % reported having bad or very bad general health
• Internal_migrant - % moved to area from England within previous year
• Dist_to_coast – nearest distance to the coast (m)
• Long – longitude
• Lat - latitude

[sebsfox]

Data loading

# get IMD from the covid tutorial
fdmr_imd <- fdmr::load_tutorial_data(
  dataset = "covid", 
  filename = "covid19_data.rds"
) |> 
  dplyr::distinct(
    MSOA11CD,
    IMD,
    perc.chinese,
    perc.indian,
    perc.bangladeshi,
    perc.pakistani,
    perc.ba,
    perc.bc,
    perc.wb
  )

# get geometry from ons geog portal
msoa11url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"
msoa11_shp <- sf::read_sf(msoa11url) |> 
  select(
    c("MSOA11CD", "geometry")
  )

df <- read.csv("data/model_data.csv") |> 
  dplyr::select(
    !c("IMD_1", starts_with("life_exp"))
  ) |> 
  dplyr::left_join(
    fdmr_imd,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  ) |> 
  dplyr::filter(
    # there are missing longs/lats (and other variables) for some areas because AreaCode is MSOA11
    !is.na(long),
    # two areas don't have an IMD
    !is.na(IMD)
  ) |> 
  left_join(
    msoa11_shp,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  )

[sebsfox]

Create mesh

df <- df |> 
  mutate(
    FYEAR = stringr::str_extract(
      FYEAR,
      "[0-9]{2}$"
    ),
    FYEAR = 2000 + as.numeric(FYEAR) - 1
  ) |> 
  sf::st_as_sf() |> 
  sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |> 
  sf::as_Spatial()

fdmr::mesh_builder(
  spatial_data = df, 
  x_coord = "long", 
  y_coord = "lat"
  )

[annabelwestermann96]

https://search.r-project.org/CRAN/refmans/fmesher/html/fm_nonconvex_hull.html

[cmmbs]

https://rpubs.com/jafet089/886687

[profrichharris]

Using a seminariogram to establish the spatial prior value

require(geostats)

tmp <- as(df, "sf")
tmp <- st_transform(tmp, 27700)
tmp$x <- st_coordinates(st_centroid(tmp))[,1]
tmp$y <- st_coordinates(st_centroid(tmp))[,2]

df2 <- tmp[tmp$FYEAR == 2014,]
semivariogram(x=df2$x,y=df2$y,z=df2$CVD)

[annabelwestermann96]

nonconvex hull mesh:

library(fmesher)

long_lat <- df |>
  select(long, lat) |>
  distinct()

sf::sf_use_s2(FALSE)
#utm  (km)
inp <- sf::st_as_sf(long_lat, coords = 1:2, crs=sf::st_crs("+proj=longlat +datum=WGS84 +no_defs"))
bnd <- fm_extensions(inp, convex = c(0.5, 2))

mesh1 <- fm_mesh_2d(boundary = bnd, max.edge = c(0.25, 2), crs=sf::st_crs("+proj=longlat +datum=WGS84 +no_defs"))

fdmr::plot_mesh(mesh1)
`

[profrichharris]

This tutorial seems quite good
https://rpubs.com/jafet089/886687

[annabelwestermann96]

library(readr)
library(dplyr)

# get IMD from the covid tutorial
fdmr_imd <- fdmr::load_tutorial_data(
  dataset = "covid", 
  filename = "covid19_data.rds"
) |> 
  dplyr::distinct(
    MSOA11CD,
    IMD,
    perc.chinese,
    perc.indian,
    perc.bangladeshi,
    perc.pakistani,
    perc.ba,
    perc.bc,
    perc.wb
  )

# get geometry from ons geog portal
msoa11url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

msoa11_shp <- sf::read_sf(msoa11url) |> 
  select(
    c("MSOA11CD", "geometry")
  )

df <- read.csv("data/model_data.csv") |> 
  dplyr::select(
    !c("IMD_1", starts_with("life_exp"))
  ) |> 
  dplyr::left_join(
    fdmr_imd,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  ) |> 
  dplyr::filter(
    # there are missing longs/lats (and other variables) for some areas because AreaCode is MSOA11
    !is.na(long),
    # two areas don't have an IMD
    !is.na(IMD)
  ) |> 
  left_join(
    msoa11_shp,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  ) |>
  mutate(
    FYEAR = stringr::str_extract(
      FYEAR,
      "[0-9]{2}$"
    ),
    FYEAR = 2000 + as.numeric(FYEAR) - 1
  )


#remotes::install_github("inlabru-org/fmesher", ref = "stable")
library(fmesher)

geo_df <- df |>
  select(long, lat, geometry, AreaCode) |>
  distinct()
  
sp_data <-  geo_df |> 
  sf::st_as_sf() |>
  sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |>
  sf::as_Spatial()

long_lat <- geo_df |>
  select(long, lat)

sf::sf_use_s2(FALSE)
#utm  (km)
inp <- sf::st_as_sf(long_lat, coords = 1:2, crs=sf::st_crs("+proj=longlat +datum=WGS84 +no_defs"))
bnd <- fm_extensions(inp, convex = c(0.5, 2))

mesh1 <- fm_mesh_2d(boundary = bnd, max.edge = c(0.25, 2), crs=sf::st_crs("+proj=longlat +datum=WGS84 +no_defs"))

fdmr::plot_mesh(mesh1)


group_index <- df$FYEAR
n_groups <- length(unique(df$FYEAR))

sp::coordinates(df) <- c("long", "lat")

# sf::df |>
#   st_as_sf() |> 
#   sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |> 
#   sf::as_Spatial()

msoa11url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

msoa11_shp <- sf::read_sf(msoa11url) |>
  sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |> 
  sf::as_Spatial()


fdmr::model_builder(spatial_data = msoa11_shp, measurement_data = df, mesh = mesh1, time_variable = "FYEAR")

[sebsfox]

For the record, I am running this from code (rather than in model builder):

initial_range <- diff(range(long_lat$long)) / 5

max_edge <- initial_range / 8


prior_range <- initial_range
spde <- INLA::inla.spde2.pcmatern(
  mesh = mesh1,
  prior.range = c(prior_range, 0.5),
  prior.sigma = c(1, 0.01)
)

rhoprior <- base::list(
  theta = list(
    prior = "pccor1", 
    param = c(0, 0.9))
)

group_index <- df$FYEAR
n_groups <- length(unique(group_index))

df@data <- df@data |> 
  mutate(
    CVD_rate = CVD / Population
  )

formula <- CVD_rate ~ 0 + Intercept(1) + IMD +
  # Resp + MSK + Psy +
  female + pop18 + pop65 +
  male +# pop_working + No_qual + Higher_qual +
  econ_inactive + bad_health + internal_migrant +
  # h_professional + 
  dist_to_coast +
  f(
    main = coordinates,
    model = spde,
    group = group_index,
    ngroup = n_groups,
    control.group = list(
      model = "ar1",
      hyper = rhoprior
    )
  )

inlabru_model <- inlabru::bru(
  formula,
  data = df,
  family = "gaussian",
  options = list(
    control.inla = list(
      reordering = "metis",
      int.strategy = "eb"
    ),
    verbose = TRUE
  )
)

[annabelwestermann96]

MSOA region lookup

https://geoportal.statistics.gov.uk/datasets/6c645fe86547424eb22783ee00bdd011_0/explore

[sebsfox]

The following code should pull together everything we've talked about, but filter for the south west region only and make bigger triangles for the mesh. The (manual) model runs in about 10 seconds - it doesn't give anything useful out, but it should be possible to iterate to a better model from there:

library(dplyr)

# get IMD from the covid tutorial
fdmr_imd <- fdmr::load_tutorial_data(
  dataset = "covid", 
  filename = "covid19_data.rds"
) |> 
  dplyr::distinct(
    MSOA11CD,
    IMD,
    perc.chinese,
    perc.indian,
    perc.bangladeshi,
    perc.pakistani,
    perc.ba,
    perc.bc,
    perc.wb
  )

# get geometry from ons geog portal
msoa11url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

msoa11_shp <- sf::read_sf(msoa11url) |> 
  select(
    c("MSOA11CD", "geometry")
  )

# msoa to region lkp
msoa_region_url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA21_BUA22_LAD22_RGN22_EW_LU_v2/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

msoa_region_lkp <- sf::read_sf(msoa_region_url) |> 
  distinct(
    MSOA21CD,
    RGN22CD,
    RGN22NM
  )

# bring in data
df <- read.csv("data/model_data.csv") |> 
  dplyr::select(
    !c("IMD_1", starts_with("life_exp"))
  ) |> 
  dplyr::left_join(
    fdmr_imd,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  ) |> 
  dplyr::filter(
    # there are missing longs/lats (and other variables) for some areas because AreaCode is MSOA11
    !is.na(long),
    # two areas don't have an IMD
    !is.na(IMD)
  ) |> 
  left_join(
    msoa11_shp,
    by = join_by(
      AreaCode == MSOA11CD
    ),
    relationship = "many-to-one"
  ) |>
  left_join(
    msoa_region_lkp,
    by = join_by(
      AreaCode == MSOA21CD
    )
  ) |>
  mutate(
    FYEAR = stringr::str_extract(
      FYEAR,
      "[0-9]{2}$"
    ),
    FYEAR = as.numeric(FYEAR),
    FYEAR = FYEAR - min(FYEAR) + 1
  ) |> 
  filter(
    RGN22NM == "South West"
  )

#remotes::install_github("inlabru-org/fmesher", ref = "stable")
library(fmesher)

geo_df <- df |>
  select(long, lat, geometry, AreaCode) |>
  distinct()

sp_data <-  geo_df |> 
  sf::st_as_sf() |>
  sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |>
  sf::as_Spatial()

long_lat <- geo_df |>
  select(long, lat)

sf::sf_use_s2(FALSE)
#utm  (km)
inp <- sf::st_as_sf(
  long_lat, 
  coords = 1:2, 
  crs=sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")
)

bnd <- fm_extensions(
  inp, 
  convex = c(0.25, 0.5)
)

mesh1 <- fm_mesh_2d(
  boundary = bnd, 
  max.edge = c(0.25, 0.5), 
  crs = sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")
)

fdmr::plot_mesh(mesh1)

group_index <- df$FYEAR
n_groups <- length(unique(df$FYEAR))

sp::coordinates(df) <- c("long", "lat")

# add CVD rate into data
df@data <- df@data |> 
  mutate(
    CVD_rate = CVD / Population
  )

msoa11url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/MSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson"

msoa11_shp <- sf::read_sf(msoa11url) |>
  filter(
    MSOA11CD %in% df$AreaCode
  ) |> 
  sf::st_transform(sf::st_crs("+proj=longlat +datum=WGS84")) |> 
  sf::as_Spatial()


# option 1: use model_builder ---------------------------------------------

fdmr::model_builder(
  spatial_data = msoa11_shp, 
  measurement_data = df, 
  mesh = mesh1, 
  time_variable = "FYEAR"
)


# option 2: manually build model -------------------------------------------------------------
initial_range <- diff(range(long_lat$long)) / 5

max_edge <- initial_range / 8


prior_range <- initial_range
spde <- INLA::inla.spde2.pcmatern(
  mesh = mesh1,
  prior.range = c(prior_range, 0.5),
  prior.sigma = c(1, 0.01)
)

rhoprior <- base::list(
  theta = list(
    prior = "pccor1", 
    param = c(0, 0.9))
)

group_index <- df$FYEAR
n_groups <- length(unique(group_index))

formula <- CVD_rate ~ 0 + Intercept(1) + IMD +
  # Resp + MSK + Psy +
  female + pop18 + pop65 +
  male +# pop_working + No_qual + Higher_qual +
  econ_inactive + bad_health + internal_migrant +
  # h_professional + 
  dist_to_coast +
  f(
    main = coordinates,
    model = spde,
    group = group_index,
    ngroup = n_groups,
    control.group = list(
      model = "ar1",
      hyper = rhoprior
    )
  )

inlabru_model <- inlabru::bru(
  formula,
  data = df,
  family = "gaussian",
  # E = df$Population,
  # control.family = list(link = "log"),
  options = list(
    control.inla = list(
      reordering = "metis",
      int.strategy = "eb"
    ),
    verbose = TRUE
  )
)

# understand model --------------------------------------------------------

fdmr::model_viewer(
  inlabru_model,
  mesh = mesh1,
  measurement_data = df
)

[annabelwestermann96]

spde <- INLA::inla.spde2.pcmatern(
                mesh = mesh,
                prior.range = c(0.05,0.1),
                prior.sigma = c(0.05,0.2)
            )

alphaprior <- list(theta = list(
                prior = 'pccor1',
                param = c(-0.2,0.8)
            ))

formula <- CVD_rate ~ 0 + Intercept(1)  +  f(
                main = coordinates,
                model = spde,
                group = group_index,
                ngroup = n_groups,
                control.group = list(
                    model = 'ar1',
                    hyper = alphaprior)
                )

model_output <- inlabru::bru(formula,
                        data = measurement_data,
                        family = 'gaussian',
                        E = NULL,
                        control.family = NULL,
                        options = list(
                            verbose = FALSE
                        )
                    )

[profrichharris]

If you want to talk to others about how you define coastal and different types of coastal...
https://www.plymouth.ac.uk/research/coastal-communities/towards-a-coastal-definition

[annabelwestermann96]

Different models to compare

spde <- INLA::inla.spde2.pcmatern(
                mesh = mesh,
                prior.range = c(0.05,0.1),
                prior.sigma = c(0.05,0.2)
            )

alphaprior <- list(theta = list(
                prior = 'pccor1',
                param = c(-0.2,0.8)
            ))

formula <- CVD_rate ~ 0 + Intercept(1)  +  f(
                main = coordinates,
                model = spde,
                group = group_index,
                ngroup = n_groups,
                control.group = list(
                    model = 'ar1',
                    hyper = alphaprior)
                )

model_output <- inlabru::bru(formula,
                        data = measurement_data,
                        family = 'gaussian',
                        E = NULL,
                        control.family = NULL,
                        options = list(
                            verbose = FALSE
                        )
                    )

[annabelwestermann96]

@sebsfox can you add yours

[profrichharris]

Before I forget, if you would like to follow-up on this work, through potentially putting together a presentation to give to some various research groups in the University (and/or work into an academic paper), or with the potential to host a PhD student for a 3-month or equivalent placement, give me shout at rich.harris@bris.ac.uk.