R4PH: Scripts

R4PH24 : 23 - 28 Sept 2024

Anthropometry

library(anthro)
library(rio)
library(eeptools)
library(ggplot2)
library(dplyr)

rm(list = ls())

## Z-scores for a single child using anthro_zscore function
anthro_zscores(sex = "M",
               age = 450,
               is_age_in_month = F,
               weight = 12,
               lenhei = 80,
               measure = "l")

## Z-scores for a data frame with many children using anthro_zscore function

# Loading the data in the environment
anthro <- rio::import("https://raw.githubusercontent.com/mgimsonline/r4ph24/refs/heads/main/data/anthro.csv") 

anthout <- anthro %>%
  mutate(age = age_calc(dob, visit_date, units = "months")) %>%
  mutate(anthro_zscores(
    sex = sex,
    age = age,
    is_age_in_month = T,
    weight = wt,
    lenhei = ht,
    measure = "h"
  ))


## Visualising the output as a density plot

anthout %>%
  ggplot() +
  geom_density(
    aes(x = zwei),
    fill = "#ff6361",
    alpha = .75,
    color = "#ff6361"
  ) +
  stat_function(data = data.frame(x = c(-4, 4)), 
                aes(x = x),
                fun = dnorm) +
  theme_minimal() +
  labs(title = "Anthropometry",
       x = "Weight-for-age z-score",
       y = element_blank())

The (not so) humble pie chart

library(dplyr)
library(ggplot2)
library(ggrepel)
library(rio)

data <- rio::import("https://raw.githubusercontent.com/mgimsonline/r4ph24/refs/heads/main/data/data.csv")

data %>% 
  group_by(phc) %>%
  count() %>%
  ungroup() %>%
  mutate(per = n/sum(n)) %>%
  ggplot(aes(x = "",
             y = per,
             fill=phc)) +
  geom_bar(position="stack", stat="identity") +
  coord_polar(theta = "y", start=0) +
  geom_label_repel(aes(label = paste0(phc," \n ",round(per*100,1),"%")),
                           position = position_stack(vjust = 0.5),
             fill = "white") +
  scale_fill_viridis_d() +
  theme_void() +
  theme(legend.position = "none") +
  labs(title = "The (not so) humble pie chart")

Heat map

library(tidyverse)
library(rio)
library(ggthemes)
library(sf)
library(rio)
library(viridis)

rm(list = ls())

wardha <- st_read(dsn = "shape_files/war_vil.shp")

data <- import("https://raw.githubusercontent.com/mgimsonline/r4ph24/refs/heads/main/data/data.csv")

cont_lims <- ggplot(data, aes(gps_long, gps_lat)) +
  geom_density_2d() +
  scale_x_continuous(limits = c(min(data[, "gps_long"]) - 2 * diff(range(data[, "gps_long"])),
                                max(data[, "gps_long"]) + 2 * diff(range(data[, "gps_long"])))) +
  scale_y_continuous(limits = c(min(data[, "gps_lat"]) - 2 * diff(range(data[, "gps_lat"])), 
                                max(data[, "gps_lat"]) + 2 * diff(range(data[, "gps_lat"]))))

cont_lims = ggplot_build(cont_lims)

rng_x <- range(cont_lims$data[[1]]$x)
x <- rng_x + c(-0.05 * diff(rng_x), 0.05 * diff(rng_x))

rng_y <- range(cont_lims$data[[1]]$y)
y <- rng_y + c(-0.05 * diff(rng_y), 0.05 * diff(rng_y))

war_talk <- wardha %>% group_by(SUB_DIST) %>% summarise()

ggplot(data) +
  geom_sf(data = wardha) +
  geom_sf(data = war_talk, aes(colour = SUB_DIST), linewidth = 1, alpha = .8) +
  geom_density_2d(aes(gps_long, gps_lat)) +
  scale_x_continuous(limits = x) +
  scale_y_continuous(limits = y) +
  stat_density_2d(data = data,
                  aes(y = gps_lat, x = gps_long, fill = ..level..),
                  geom = "polygon") +
  scale_fill_viridis(alpha = .3, direction = -1) +
  coord_sf(xlim = c(78.3, 78.9), ylim = c(20.5, 21.0)) +
  theme_map() +
  theme(legend.position = "none") +
  labs(title = "Population density of Sevagram HDSS")

Loops

library(rio)
library(dplyr)
library(gt)

data <- import("https://raw.githubusercontent.com/mgimsonline/r4ph24/refs/heads/main/data/data.csv")

## For loops
for (i in unique(data$phc)) {
  plot <- data %>%
    filter(phc == i) %>%
    group_by(stress) %>%
    summarise(sbp = mean(sbp, na.rm = T)) %>%
    gt() %>%
    tab_header(title = paste0("Sys BP by stress for ", i, " PHC")) %>%
    gtsave(paste0(i, "_stress_sbp.pdf"))
}

## While loops
x = 1
while (x < 10000000000000000) {
  x = x*2
  print(as.character(x))
}

## If else
if(mean(data$sbp) > 120) {
  print("Mean SBP is more than 120")
} else {
  print("Mean SBP is less than 120")
}

Kobotoolbox

library(KoboconnectR)
library(xml2)
library(XML)
library(dplyr)
library(gt)

rm(list = ls())

url <- "kf.kobotoolbox.org" ## Change the server URL
username <- "username" ## Change the user name
password <- "password" ## Change the password
project <- "Project Name" ## Change the project name

## Get a list of projects
assets <- kobotools_api(url = url,
                        uname = username,
                        pwd = password)

## Select the project
id <- assets %>% 
  filter(name == project) %>% 
  .$assetid

## Request the data from the server
data <- kobotools_kpi_data(
  url = url,
  uname = username,
  pwd = password,
  assetid = id,
  format = "xml"
)

## Convert the data to data frame
xml <- xmlParse(data)

root_node <- xmlRoot(xml)

xml_results <- getNodeSet(root_node, "//results")[[1]]

results <- xmlToDataFrame(xml_results)