NDMA_Deidentification_Jul2024
Pheroze
2024-10-18
NDMA Data De-Identification
The following are steps undertaken for deidentifying NDMA data. This dataset covers the month of July 2024 for HHA, KIA and MUAC. Since the geocoordinates in the HHA dataset represent household coordinates, we will mask them (random displacement) using the Haversine Formula to randomly distribute a point around a central coordinate within a radius of 2.5 KM and drop other P.I.I.s in the relevant sheets.
1. HHA
# Load required libraries
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.2library(geosphere)## Warning: package 'geosphere' was built under R version 4.3.3## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)library(openxlsx)## Warning: package 'openxlsx' was built under R version 4.3.3file_path <- "C:/Users/AAH USER/Downloads/HHA, KIA & MUAC July 2024.xlsx"
# Read the specific HHA sheet into a data frame
hha_jul_2024 <- read.xlsx(file_path, sheet = "HHA", startRow = 2)The dataset contains household coordinates in columns “Lat” and “Long” which are considered P.I.I’s so we mask the coordinates, verify by plotting a histogram of the distribution of displacement distances of the original and displaced coordinates to establish uniformity.
We check for and deal with outliers if any in the “Lat” and “Long” columns
# Verify the dataset
summary(hha_jul_2024$Lat)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -4.45349 -1.70092 -0.02985 -0.08483 1.41839 22.32788 30summary(hha_jul_2024$Long)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34.88 36.78 37.88 38.02 39.50 91.83 30There seems to be an erroneous entry since the max values in the “Lat” and “Long” columns are well above the mean. 30 rows also dont have entries for coordinates. We deal with this entry.
# Calculate the mean for Lat and Long, ignoring NA values
lat_mean <- mean(hha_jul_2024$Lat, na.rm = TRUE)
lon_mean <- mean(hha_jul_2024$Long, na.rm = TRUE)
lat_sd <- sd(hha_jul_2024$Lat, na.rm = TRUE)
lon_sd <- sd(hha_jul_2024$Long, na.rm = TRUE)
# Calculate Z-scores
hha_jul_2024 <- hha_jul_2024 %>%
mutate(lat_z = (Lat - lat_mean) / lat_sd,
lon_z = (Long - lon_mean) / lon_sd)
# Set threshold for identifying outliers
threshold <- 3 # Common threshold for Z-scores
# Replace outliers with NA
hha_jul_2024 <- hha_jul_2024 %>%
mutate(Lat = ifelse(abs(lat_z) > threshold & !is.na(lat_z), NA, Lat),
Long = ifelse(abs(lon_z) > threshold & !is.na(lon_z), NA, Long))
# Remove the Z-score columns
hha_jul_2024 <- hha_jul_2024 %>%
select(-lat_z, -lon_z)
# Verify the dataset
summary(hha_jul_2024$Lat)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -4.45349 -1.70092 -0.02986 -0.09459 1.41837 4.08103 32summary(hha_jul_2024$Long)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34.88 36.78 37.87 38.00 39.50 41.86 32Only two rows are affected by the outlier and are replaced with NAs, bringing the total to 32 NAs. We proceed to mask the coordinates.
# Create backup columns for original coordinates
hha_jul_2024$Original_Lat <- hha_jul_2024$Lat
hha_jul_2024$Original_Long <- hha_jul_2024$Long
# Function to generate random displaced coordinates with uniform distance distribution
mask_coordinates_uniform <- function(lat, lon, radius_km) {
R <- 6371 # Earth radius in kilometers
# Random bearing angle (in radians)
bearing <- runif(1, 0, 2 * pi)
# Random distance uniformly sampled from [0, radius_km]
rand_dist <- runif(1, 0, radius_km) / R # Uniformly sampled distance in radians
# Convert original coordinates to radians
lat_rad <- lat * pi / 180
lon_rad <- lon * pi / 180
# Calculate new latitude (in radians)
new_lat <- asin(sin(lat_rad) * cos(rand_dist) +
cos(lat_rad) * sin(rand_dist) * cos(bearing))
# Calculate new longitude (in radians)
new_lon <- lon_rad + atan2(sin(bearing) * sin(rand_dist) * cos(lat_rad),
cos(rand_dist) - sin(lat_rad) * sin(new_lat))
# Convert back to degrees
new_lat <- new_lat * 180 / pi
new_lon <- new_lon * 180 / pi
return(c(new_lat, new_lon))
}
# Set displacement radius in kilometers
radius_km <- 2.5
# Generate masked coordinates for each row using the modified function
masked_coords <- t(apply(hha_jul_2024, 1, function(row) {
mask_coordinates_uniform(as.numeric(row["Original_Lat"]), as.numeric(row["Original_Long"]), radius_km)
}))
# Replace the original Lat and Long columns with the masked coordinates
hha_jul_2024$Lat <- masked_coords[, 1]
hha_jul_2024$Long <- masked_coords[, 2]We then evaluate the distribution of the displacement distances before dropping the original coordinates by plotting a histogram and a scatter plot.
# Calculate displacement distances (in kilometers) as before
displacement_distances <- distHaversine(
cbind(hha_jul_2024$Long, hha_jul_2024$Lat), # Masked coordinates
cbind(hha_jul_2024$Original_Long, hha_jul_2024$Original_Lat) # Original coordinates
) / 1000 # Convert meters to kilometers
# Add displacement distances to the dataset for further analysis
hha_jul_2024$Displacement_Distance <- displacement_distances
# Plot: Histogram of displacement distances
ggplot(hha_jul_2024, aes(x = Displacement_Distance)) +
geom_histogram(binwidth = 0.1, fill = "skyblue", color = "black") +
labs(title = "Distribution of Displacement Distances",
x = "Displacement Distance (km)", y = "Count") +
theme_minimal()## Warning: Removed 32 rows containing non-finite values (`stat_bin()`).
There are no significant peaks or valleys in the histogram, suggesting that the displacements are indeed more uniformly distributed, as intended.
# Plot: Original vs Masked Coordinates Scatterplot
ggplot() +
geom_point(data = hha_jul_2024, aes(x = Original_Long, y = Original_Lat, color = "Original"),
alpha = 0.5, size = 1, shape = 16) +
geom_point(data = hha_jul_2024, aes(x = Long, y = Lat, color = "Displaced"),
alpha = 0.5, size = 1, shape = 16) +
labs(title = "Original vs Displaced Coordinates",
x = "Longitude", y = "Latitude", color = "Coordinate Type") + # Legend title
scale_color_manual(values = c("Original" = "red", "Displaced" = "blue")) + # Custom colors
theme_minimal() +
coord_fixed() # Ensures aspect ratio is 1:1 for accurate geographic representation## Warning: Removed 32 rows containing missing values (`geom_point()`).
## Removed 32 rows containing missing values (`geom_point()`).
The scatterplot displays the original coordinates in red and the displaced coordinates in blue. It helps visualize how the coordinates were shifted.
We then drop the original coordinates column leaving only the masked coordinates columns. We also drop other PII’s which are the “HouseholdName”, “HouseHoldHead”, and “RespondentName”.
# Drop the specified PII columns along with original coordinates
hha_jul_2024 <- hha_jul_2024 %>%
select(-c(Original_Lat, Original_Long, Displacement_Distance, HouseholdName, HouseHoldHead, RespondentName))
# Check the updated dataset
head(hha_jul_2024)## QID County SubCounty Ward LivelihoodZone Month Year Lat Long
## 1 519096 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2738843 37.55568
## 2 519097 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2686389 37.56314
## 3 519098 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2888230 37.53497
## 4 519099 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2587653 37.56396
## 5 519101 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2680784 37.55187
## 6 519102 Isiolo Isiolo North Burat Agro Pastoral July 2024 0.2716081 37.55187
## InterviewDate HouseholdCode HeadEducationLevel
## 1 45481.00 ISL0421 1. Primary
## 2 45482.00 ISL0430 1. Primary
## 3 45482.82 ISL0419 1. Primary
## 4 45483.00 ISL0428 None
## 5 45483.00 ISL0427 None
## 6 45483.37 ISL0426 None
## MainHHIncomeSource HeadGender RespondentGender
## 1 2. Sale of livestock and livestock products Female Female
## 2 3. Sale of crops Female Female
## 3 4. Casual labour Female Female
## 4 2. Sale of livestock and livestock products Female Female
## 5 2. Sale of livestock and livestock products Female Female
## 6 2. Sale of livestock and livestock products Female Female
## MaleMembers FemaleMembers ChildrenBelow5 KeepLivestock MilkAnimals
## 1 4 4 2 TRUE TRUE
## 2 4 3 3 TRUE FALSE
## 3 3 2 2 TRUE FALSE
## 4 3 3 2 TRUE FALSE
## 5 3 2 3 TRUE FALSE
## 6 5 3 2 TRUE FALSE
## MilkSource HowOftenMilked AverageMilkedPerDay
## 1 Own livestock production NA NA
## 2 NA NA
## 3 <NA> NA NA
## 4 NA NA
## 5 NA NA
## 6 <NA> NA NA
## AverageMilkConsumedPerDay WhoDrankMilk AverageMilkPrice
## 1 NA Children under 5 years NA
## 2 NA NA
## 3 NA <NA> NA
## 4 NA Children under 5 years NA
## 5 NA Children under 5 years NA
## 6 NA Children under 5 years NA
## HarvestedInLastWeeks AcresHarvested BagsHarvested HaveFoodStock
## 1 FALSE NA NA TRUE
## 2 FALSE NA NA FALSE
## 3 TRUE NA NA FALSE
## 4 FALSE NA NA TRUE
## 5 FALSE NA NA FALSE
## 6 TRUE NA NA FALSE
## FoodStockSources DaysStockLast WaterSource1 WaterSource2 WaterSource3
## 1 Purchase 3 Rivers
## 2 NA Rivers
## 3 <NA> NA Rivers <NA> <NA>
## 4 Purchase 3 Rivers
## 5 NA Rivers
## 6 <NA> NA Rivers <NA> <NA>
## NormalWaterSource WhyNotNormalWaterSource DaysWaterSourceExpectedToLast
## 1 TRUE Breakdown of water source 12
## 2 TRUE Breakdown of water source 12
## 3 TRUE <NA> 12
## 4 TRUE Breakdown of water source 12
## 5 TRUE Breakdown of water source 12
## 6 TRUE <NA> 12
## DistanceFromWaterSource NoWaterJerryCans JerryCansCost
## 1 2.6 10 NA
## 2 2.5 9 NA
## 3 1.5 8 NA
## 4 2.5 8 NA
## 5 2.5 8 NA
## 6 2.6 10 NA
## NormalHHWaterConsumption HHPayForWater CostTransportJerryCan
## 1 6 FALSE NA
## 2 5 FALSE NA
## 3 5 FALSE NA
## 4 5 FALSE NA
## 5 5 FALSE NA
## 6 5 FALSE NA
## TreatWaterBeforeDrinking WaterTreatmentMethodUsed CSI_ReliedOnLess
## 1 FALSE 2
## 2 FALSE 2
## 3 FALSE <NA> 2
## 4 FALSE 2
## 5 FALSE 3
## 6 FALSE <NA> 1
## CSI_BorrowedFood CSI_ReducedNoOfMeals CSI_ReducedPortionMealSize
## 1 0 1 0
## 2 0 3 1
## 3 0 3 2
## 4 0 2 1
## 5 0 3 1
## 6 0 2 2
## CSI_QuantityForAdult CSI_SoldHouseholdAssets CSI_ReducedNonFoodExpenses
## 1 0 1 1
## 2 0 1 1
## 3 2 1 1
## 4 0 1 1
## 5 0 1 1
## 6 3 1 1
## CSI_SoldProductiveAssets CSI_SpentSavings CSI_BorrowedMoney CSI_SoldHouseLand
## 1 1 1 1 1
## 2 1 1 1 1
## 3 1 1 1 1
## 4 1 1 1 1
## 5 1 1 1 1
## 6 1 1 1 1
## CSI_WithdrewChildrenSchool CSI_SoldLastFemaleAnimal CSI_Begging
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1
## 5 1 1 1
## 6 1 1 1
## CSI_SoldMoreAnimals HFC_GrainDays HFC_GrainSource HFC_RootsDays
## 1 1 7 5 NA
## 2 1 7 5 NA
## 3 1 7 5 NA
## 4 1 7 5 NA
## 5 1 7 5 NA
## 6 1 7 5 NA
## HFC_RootsSource HFC_PulsesNutsDays HFC_PulsesNutsSource HFC_OrangeVegDays
## 1 NA 4 5 NA
## 2 NA 5 5 NA
## 3 NA 0 NA NA
## 4 NA 3 5 NA
## 5 NA 3 5 NA
## 6 NA 0 NA NA
## HFC_OrangeVegSource HFC_GreenLeafyDays HFC_GreenLeafySource HFC_OtherVegDays
## 1 NA NA NA 5
## 2 NA NA NA 4
## 3 NA NA NA 2
## 4 NA NA NA 5
## 5 NA NA NA 2
## 6 NA NA NA 3
## HFC_OtherVegSource HFC_OrangeFruitsDays HFC_OrangeFruitsSource
## 1 5 NA NA
## 2 5 NA NA
## 3 5 NA NA
## 4 1 NA NA
## 5 5 NA NA
## 6 5 NA NA
## HFC_OtherFruitsDays HFC_OtherFruitsSource HFC_MeatDays HFC_MeatSource
## 1 3 1 0 10
## 2 3 5 0 10
## 3 3 5 0 NA
## 4 3 5 0 10
## 5 1 5 0 10
## 6 2 5 0 NA
## HFC_LiverDays HFC_LiverSource HFC_FishDays HFC_EggsDays HFC_EggsSource
## 1 NA NA NA NA NA
## 2 NA NA NA NA NA
## 3 NA NA NA NA NA
## 4 NA NA NA NA NA
## 5 NA NA NA NA NA
## 6 NA NA NA NA NA
## HFC_MilkDays HFC_MilkSource HFC_OilDays HFC_OilSource HFC_SugarDays
## 1 2 5 6 5 6
## 2 2 1 7 5 7
## 3 2 5 6 5 7
## 4 2 1 7 5 7
## 5 0 10 7 5 7
## 6 2 1 7 5 6
## HFC_SugarSource HFC_CondimentsDays HFC_CondimentsSource MainIncomeSource
## 1 5 0 NA 1. Sale of crops
## 2 5 0 NA 1. Sale of crops
## 3 5 0 NA 4. Casual labour
## 4 5 3 5 5. Sale of wood
## 5 5 3 5 4. Casual labour
## 6 5 0 NA 6. Sale of charcoal
## MaleCasualLabour FemaleCasualLabour CasualLabourEarn CharcoalSaleEarn
## 1 1 1 2000 NA
## 2 NA NA NA NA
## 3 1 NA 1200 NA
## 4 NA NA NA NA
## 5 1 1 1200 NA
## 6 NA NA NA NA
## WoodSaleEarn DivisionID CountyID SiteID LivelihoodZoneID DateCaptured
## 1 NA 39 4 1 2 NA
## 2 NA 39 4 1 2 NA
## 3 NA 39 4 1 2 NA
## 4 NA 39 4 1 2 NA
## 5 NA 39 4 1 2 NA
## 6 NA 39 4 1 2 NAWe also have to ensure that the “InterviewDate” column is parsed correctly as a date before saving the worksheet to the new workbook.
# Ensure the column is numeric
hha_jul_2024$InterviewDate <- as.numeric(hha_jul_2024$InterviewDate)
# Convert the numeric date to Date format
hha_jul_2024$InterviewDate <- as.Date(hha_jul_2024$InterviewDate, origin = "1899-12-30")
# View the first few dates to verify the conversion
head(hha_jul_2024$InterviewDate)## [1] "2024-07-08" "2024-07-09" "2024-07-09" "2024-07-10" "2024-07-10"
## [6] "2024-07-10"We save the sheet to a new workbook for the deidentified data for NDMA July 2024. We will write subsequent deidentified KIA and MUAC sheets to this workbook.
# Define the path for the new workbook
new_file_path <- "C:/Users/AAH USER/OneDrive - Action Against Hunger USA/Documents/NDMA_DeIdentified/HHA_KIA_MUAC_July_2024.xlsx"
# Create a new workbook
wb <- createWorkbook()
# Add the HHA data to the new workbook
addWorksheet(wb, "HHA")
writeData(wb, "HHA", hha_jul_2024)
# Save the new workbook
saveWorkbook(wb, new_file_path, overwrite = TRUE)2. KIA
# Load required libraries
library(dplyr)
library(ggplot2)
library(geosphere)
library(openxlsx)
file_path <- "C:/Users/AAH USER/Downloads/HHA, KIA & MUAC July 2024.xlsx"
# Read the specific HHA sheet into a data frame
kia_jul_2024 <- read.xlsx(file_path, sheet = "KIA", startRow = 2)The dataset contains coordinates, however these are not household coordinates so we do not need to mask them. We ensure the Interview date column is parsed correctly and save this subsequent sheet to the workbook.
# Ensure the column is numeric
kia_jul_2024$InterviewDate <- as.numeric(kia_jul_2024$InterviewDate)
# Convert the numeric date to Date format
kia_jul_2024$InterviewDate <- as.Date(kia_jul_2024$InterviewDate, origin = "1899-12-30")
# View the first few dates to verify the conversion
head(kia_jul_2024$InterviewDate)## [1] "2024-07-12" "2024-07-12" "2024-07-11" "2024-07-12" "2024-07-11"
## [6] "2024-07-13"We save this new sheet alongside the previous in the workbook created
# Define the path for the existing Excel workbook
existing_file_path <- "C:/Users/AAH USER/OneDrive - Action Against Hunger USA/Documents/NDMA_DeIdentified/HHA_KIA_MUAC_July_2024.xlsx"
# Load the existing workbook
wb <- loadWorkbook(existing_file_path)
# Add the KIA data to the existing workbook
addWorksheet(wb, "KIA ")
writeData(wb, "KIA ", kia_jul_2024)
# Save the updated workbook
saveWorkbook(wb, existing_file_path, overwrite = TRUE)3. MUAC
# Load required libraries
library(dplyr)
library(ggplot2)
library(geosphere)
library(openxlsx)
file_path <- "C:/Users/AAH USER/Downloads/HHA, KIA & MUAC July 2024.xlsx"
# Read the specific HHA sheet into a data frame
muac_jul_2024 <- read.xlsx(file_path, sheet = "MUAC", startRow = 2)This data set has P.I.I’s in the “ChildName” column so we will drop that, ensure the Interview date column is parsed correctly and save this subsequent sheet to the workbook.
# Drop the specified PII column
muac_jul_2024 <- muac_jul_2024 %>%
select(-ChildName)
# Ensure the column is numeric
muac_jul_2024$InterviewDate <- as.numeric(muac_jul_2024$InterviewDate)
# Convert the numeric date to Date format
muac_jul_2024$InterviewDate <- as.Date(muac_jul_2024$InterviewDate, origin = "1899-12-30")
# Check the updated dataset
head(muac_jul_2024)## MUACIndicatorID QID County SubCounty Ward LivelihoodZone Month
## 1 1440342 515609 Kilifi Kilifi South Mwarakaya Food Cropping July
## 2 1440343 515609 Kilifi Kilifi South Mwarakaya Food Cropping July
## 3 1440344 515609 Kilifi Kilifi South Mwarakaya Food Cropping July
## 4 1440345 515609 Kilifi Kilifi South Mwarakaya Food Cropping July
## 5 1440346 515611 Kilifi Kilifi South Mwarakaya Food Cropping July
## 6 1440347 515611 Kilifi Kilifi South Mwarakaya Food Cropping July
## Year HouseholdCode Gender MUAC MUAC_Color AgeInMonths LiveInHousehold
## 1 2024 KIL0703 Female 153 NA 27 TRUE
## 2 2024 KIL0703 Female 152 NA 40 FALSE
## 3 2024 KIL0703 Male 152 NA 34 FALSE
## 4 2024 KIL0703 Male 152 NA 24 FALSE
## 5 2024 KIL0701 Male 151 NA 57 TRUE
## 6 2024 KIL0701 Male 150 NA 59 FALSE
## SufferedIllnesses InterviewDate DivisionID CountyID SiteID LivelihoodZoneID
## 1 <NA> 2024-07-01 90 15 189 1019
## 2 <NA> 2024-07-01 90 15 189 1019
## 3 <NA> 2024-07-01 90 15 189 1019
## 4 <NA> 2024-07-01 90 15 189 1019
## 5 <NA> 2024-07-01 90 15 189 1019
## 6 <NA> 2024-07-01 90 15 189 1019We save this new sheet alongside the previous two in the workbook created
# Define the path for the existing Excel workbook
existing_file_path <- "C:/Users/AAH USER/OneDrive - Action Against Hunger USA/Documents/NDMA_DeIdentified/HHA_KIA_MUAC_July_2024.xlsx"
# Load the existing workbook
wb <- loadWorkbook(existing_file_path)
# Add the MUAC data to the existing workbook
addWorksheet(wb, "MUAC ")
writeData(wb, "MUAC ", muac_jul_2024)
# Save the updated workbook
saveWorkbook(wb, existing_file_path, overwrite = TRUE)