We chose to model the deployment of Method X, a contraceptive implant, in Kenya. For the purpose of interpreting results later on, we assume that Method X offers a unique balance between short-term and long-term contraceptive options. Unlike daily pills that require strict adherence, Method X provides long-lasting protection while avoiding the extended commitment of traditional implants (which typically last three years) or IUDs (which can last five or more years).

Additionally, Method X is minimally invasive, requiring only a simple clinical insertion rather than a surgical procedure. Unlike pills, which rely on daily user compliance, Method X eliminates the risk of missed doses, making it a more reliable option. At the same time, it does not require the same level of permanence as other long-term methods, allowing women greater flexibility in their reproductive choices.

Since Method X requires administration at a clinic, access to healthcare services plays a crucial role in its adoption. This makes it particularly relevant for evaluating disparities in contraceptive access between rural and urban populations. By modeling adoption trends, we can better understand the factors influencing its uptake, including income level, marital status, and prior contraceptive use.

  1. Loading Necessary Libraries

We use library dplyr to help in handling and processing data, such as filtering, summarizing, and modifying simulated data. The library ggplot2 is essential for creating clear and professional visualizations, like trend lines for adoption rates over time. We also considered tidyr which allows us to reshape and organize data efficiently, such as pivoting wide-format data into long format for plotting multiple scenarios. Together, these libraries streamline data analysis and make it easier to interpret and present your results.

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, union
library(ggplot2)
library(tidyr)

Simulating Contraceptive Use and Fertility Data for Women in Kenya (Ages 15-49)

This simulation models 100,000 women aged 15 to 49, focusing on reproductive health factors such as contraceptive use, marital status, number of children, and access to family planning services. The probabilities assigned to different variables are based on real-world data trends from sources like the Kenya Demographic and Health Survey (KDHS) 2022.

1. Age Distribution

The dataset reflects the typical population structure in Kenya, where younger women (15–24) make up a larger proportion, while the number of women gradually declines in older age groups (35–49).

2. Number of Children

The likelihood of having children increases with age. Teenagers (15–19) are least likely to have children, while women aged 30 and above generally have multiple children. The average number of children in the dataset is approximately 3.4, reflecting common fertility patterns.

3. Marital Status

Marital status is strongly influenced by age:
- Women under 20 are almost always single.
- Between 20–24, marriage becomes slightly more likely but remains uncommon.
- From 25–29, marriage rates increase significantly.
- By 30–34, most women are married.
- From 35–49, marriage is at its highest prevalence.

This aligns with the reality that women in Kenya typically marry later than in previous generations, influenced by education, career growth, and economic factors.

4. Location (Rural vs. Urban)

Younger women are more likely to live in urban areas, while older women are more likely to be in rural settings. This reflects trends where younger populations migrate to cities for work and education, while older women settle in more rural environments.

5. Income Levels

Income level is distributed based on age and location:
- Younger women and those in rural areas tend to be in the low-income bracket.
- Middle-aged women (30–39) see a shift towards middle-income levels.
- Older women (40–49), especially those in urban areas, have higher incomes.

6. Access to Family Planning Services

Access to clinics determines whether a woman can obtain contraceptive implants (Method X).
- Women with higher incomes have the greatest access.
- Women in middle-income brackets have moderate access.
- Low-income and rural women face limited access to family planning services.

This reflects real-world disparities where economic status and location affect healthcare access.

7. Contraceptive Use & Campaign Impact

The dataset establishes a baseline where most women initially use “Other” methods (e.g., pills, natural methods, or no contraception). A campaign is introduced to promote Method X (an implant requiring clinical insertion), and the progression of adoption is tracked.

Accidental pregnancies are also included, as they may influence a woman’s decision to switch to Method X after experiencing an unintended pregnancy.

Conclusion

This simulation provides a realistic representation of contraceptive use and reproductive health among Kenyan women. By modeling key factors such as age, marital status, income, and access to clinics, it helps analyze how different groups respond to family planning campaigns and which women are more likely to adopt Method X over time.

set.seed(108)
ages <- 15:49  
prob1 <- c(
  rep(0.22 / 5, 5),  # 15-19
  rep(0.19 / 5, 5),  # 20-24
  rep(0.17 / 5, 5),  # 25-29
  rep(0.15 / 5, 5),  # 30-34
  rep(0.10 / 5, 5),  # 35-39
  rep(0.095 / 5, 5),  # 40-44
  rep(0.075 / 5, 5)  # 45-49
)
prob2 <- c(1,0)

sampled_ages <- sample(ages, 100000, replace = TRUE, prob = prob1)
cmethods <- c("Others","Method X")

# Define number of children probabilities based on age group
children_probs <- list(
  "15-19" = c(0.85, 0.15, 0, 0, 0, 0),  
  "20-24" = c(0.75, 0.25, 0, 0, 0, 0),  
  "25-29" = c(0.3, 0.45, 0.15, 0, 0, 0),  
  "30-34" = c(0.02, 0.25, 0.45, 0.2, 0.07, 0.01),  
  "35-39" = c(0.01, 0.02, 0.3, 0.4, 0.2, 0.07),  
  "40-44" = c(0.005, 0.01, 0.1, 0.35, 0.35, 0.085),  
  "45-49" = c(0.002, 0.005, 0.02, 0.1, 0.38, 0.493)  
)

assign_children <- function(age) {
  if (age < 20) {
    return(sample(1:6, 1, prob = children_probs[["15-19"]]))
  } else if (age >= 20 & age < 25) {
    return(sample(1:6, 1, prob = children_probs[["20-24"]]))
  } else if (age >= 25 & age < 30) {
    return(sample(1:6, 1, prob = children_probs[["25-29"]]))
  } else if (age >= 30 & age < 35) {
    return(sample(1:6, 1, prob = children_probs[["30-34"]]))
  } else if (age >= 35 & age < 40) {
    return(sample(1:6, 1, prob = children_probs[["35-39"]]))
  } else if (age >= 40 & age < 45) {
    return(sample(1:6, 1, prob = children_probs[["40-44"]]))
  } else {
    return(sample(1:6, 1, prob = children_probs[["45-49"]]))
  }
}

num_children <- sapply(sampled_ages, assign_children)

marital_status <- sapply(sampled_ages, function(age) {
  if (age < 20) {
    return("Single")
  } else if (age >= 20 & age < 25) {
    return(sample(c("Married", "Single"), 1, prob = c(0.2, 0.8)))
  } else if (age >= 25 & age < 30) {
    return(sample(c("Married", "Single"), 1, prob = c(0.5, 0.5)))
  } else if (age >= 30 & age < 35) {
    return(sample(c("Married", "Single"), 1, prob = c(0.75, 0.25)))
  } else {
    return(sample(c("Married", "Single"), 1, prob = c(0.9, 0.1)))
  }
})

location <- ifelse(sampled_ages < 30, 
                   sample(c("Urban", "Rural"), 100000, replace = TRUE, prob = c(0.63, 0.36)), 
                   sample(c("Urban", "Rural"), 100000, replace = TRUE, prob = c(0.73, 0.27)))

income_level <- sapply(sampled_ages, function(age) {
  if (age < 20) {
    return(sample(c("Low","Middle","High"), 1, prob = c(0.69, 0.26,0.05)))
  } else if (age >= 20 & age < 25) {
    return(sample(c("Low","Middle","High"), 1, prob = c(0.77, 0.21,0.02)))
  } else if (age >= 25 & age < 30) {
    return(sample(c("Low","Middle","High"), 1, prob = c(0.53, 0.31,0.16)))
  } else if (age >= 30 & age < 35) {
    return(sample(c("Low","Middle","High"), 1, prob = c(0.27, 0.61,0.12)))
  } else {
    return(sample(c("Low","Middle","High"), 1, prob = c(0.21, 0.63,0.16)))
  }
})

access_family_planning <- ifelse(income_level == "High", "Yes", 
                          ifelse(income_level == "Middle", sample(c("Yes", "No"), 100000, replace = TRUE, prob = c(0.89, 0.11)),
                                 sample(c("Yes", "No"), 100000, replace = TRUE, prob = c(0.74, 0.26))))

# Assign education level based on age and income
education_level <- sapply(1:100000, function(i) {
  age <- sampled_ages[i]
  income <- income_level[i]
  
  if (age < 20) {
    return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.55, 0.43, 0.02)))
  } else if (age >= 20 & age < 25) {
    if (income == "Low") {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.50, 0.45, 0.05)))
    } else {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.35, 0.50, 0.15)))
    }
  } else if (age >= 25 & age < 30) {
    if (income == "Low") {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.40, 0.50, 0.10)))
    } else {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.25, 0.50, 0.25)))
    }
  } else if (age >= 30 & age < 35) {
    if (income == "Low") {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.30, 0.50, 0.20)))
    } else {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.15, 0.50, 0.35)))
    }
  } else {
    if (income == "Low") {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.25, 0.50, 0.25)))
    } else {
      return(sample(c("Primary", "Secondary", "Higher"), 1, prob = c(0.10, 0.50, 0.40)))
    }
  }
})

# Create the final population data frame
population <- data.frame(
  id = 1:100000,
  age = sampled_ages,
  current_method = sample(cmethods, 100000, replace = TRUE, prob = prob2),
  num_children = num_children,
  marital_status = marital_status,
  location = location,
  income_level = income_level,
  access_family_planning = access_family_planning,
  education_level = education_level
)


#View(population)

head(population)
##   id age current_method num_children marital_status location income_level
## 1  1  27         Others            1         Single    Rural       Middle
## 2  2  24         Others            2         Single    Urban          Low
## 3  3  20         Others            2         Single    Urban          Low
## 4  4  32         Others            4        Married    Rural       Middle
## 5  5  27         Others            2        Married    Rural         High
## 6  6  17         Others            1         Single    Rural       Middle
##   access_family_planning education_level
## 1                    Yes       Secondary
## 2                    Yes       Secondary
## 3                    Yes       Secondary
## 4                    Yes         Primary
## 5                    Yes       Secondary
## 6                    Yes         Primary
#install.packages("NetLogoR")
library(NetLogoR)
library(dplyr)
library(NetLogoR)

# Create a properly initialized world
world <- createWorld(minPxcor = 0, maxPxcor = 100, minPycor = 0, maxPycor = 100)

# Check world properties
print(world)
## class       : worldMatrix 
## resolution  : 1, 1 (x, y)
## dimensions  : Pxcor:  0 , 100 
##               Pycor:  0 , 100 
## names       : layer 
## min values  :    NA 
## max values  :    NA 
## First 4 rows and  4 columns:
##      [,1] [,2] [,3] [,4]
## [1,]   NA   NA   NA   NA
## [2,]   NA   NA   NA   NA
## [3,]   NA   NA   NA   NA
## [4,]   NA   NA   NA   NA
print(maxPxcor(world))  
## [1] 100
# Number of agents
num_agents <- 100000
sampled_ages <- sample(35:49, num_agents, replace = TRUE)
cmethods <- c("Others", "Implant", "None")  # Example contraceptive methods
prob2 <- c(0.5, 0.3, 0.2)  # Example probabilities for methods

num_children <- sample(0:5, num_agents, replace = TRUE)  # Example children count
marital_status <- sample(c("Single", "Married"), num_agents, replace = TRUE)
location <- sample(c("Urban", "Rural"), num_agents, replace = TRUE)
income_level <- sample(c("Low", "Medium", "High"), num_agents, replace = TRUE)
access_family_planning <- sample(c(TRUE, FALSE), num_agents, replace = TRUE)
education_level <- sample(c("None", "Primary", "Secondary", "Higher"), num_agents, replace = TRUE)
# Creating agent population data frame
population <- data.frame(
  id = 1:num_agents,
  age = sampled_ages,
  current_method = sample(cmethods, num_agents, replace = TRUE, prob = prob2),
  num_children = num_children,
  marital_status = marital_status,
  location = location,
  income_level = income_level,
  access_family_planning = access_family_planning,
  education_level = education_level
)
View(population)
# Create agents (turtles) in the world
agents <- createTurtles(n = num_agents, world = world)

# Assign attributes directly to agentMatrix
agents@.Data <- cbind(
  agents@.Data,
  age = population$age,
  current_method = population$current_method,
  num_children = population$num_children,
  marital_status = population$marital_status,
  location = population$location,
  income_level = population$income_level,
  access_family_planning = population$access_family_planning,
  education_level = population$education_level
)
agents@.Data <- agents@.Data[, !(colnames(agents@.Data) %in% c("prevX", "prevY", "breed", "color","xcor","ycor","heading","who"))]
# Check agents
#print(agents@.Data)
#str(agents@.Data)  # Check structure of agents
head(agents@.Data)  # View first few rows
##      age  current_method num_children marital_status location income_level
## [1,] "39" "None"         "5"          "Married"      "Urban"  "High"      
## [2,] "44" "None"         "5"          "Single"       "Urban"  "Low"       
## [3,] "46" "Implant"      "3"          "Married"      "Rural"  "High"      
## [4,] "45" "Others"       "3"          "Single"       "Rural"  "Medium"    
## [5,] "48" "Others"       "0"          "Single"       "Urban"  "Low"       
## [6,] "35" "None"         "1"          "Single"       "Rural"  "High"      
##      access_family_planning education_level
## [1,] "TRUE"                 "Primary"      
## [2,] "TRUE"                 "Higher"       
## [3,] "FALSE"                "Primary"      
## [4,] "FALSE"                "Secondary"    
## [5,] "FALSE"                "Secondary"    
## [6,] "FALSE"                "Higher"
table(agents@.Data[, "current_method"])  # Check distribution of current_method
## 
## Implant    None  Others 
##   29763   19984   50253
# Define adoption function
adopt_method_x <- function(agents) {
  adoption_rate <- sample(seq(0.03, 0.07, by=0.01), 1)  # Random rate between 3-7%
  
  # Convert current_method column to character to avoid factor issues
  current_methods <- agents@.Data[, "current_method"]
  
  # Find indices of agents using "None" (not "Others")
  none_indices <- which(current_methods == "None")
  
  # Select subset that adopts Method X
  adopt_indices <- none_indices[runif(length(none_indices)) < adoption_rate]
  
  # Apply changes in-place
  agents@.Data[adopt_indices, "current_method"] <- "Method X"
  
  return(agents)
}


# Define discontinuation function
discontinue_method_x <- function(agents, high_side_effects = FALSE) {
  base_rate <- 0.05  # 5% default discontinuation
  if (high_side_effects) base_rate <- base_rate * 2  # Double if high side effects
  
  # Find agents currently using "Method X"
  method_x_idx <- which(agents@.Data[, "current_method"] == "Method X")
  
  # Randomly determine who discontinues
  discontinue_idx <- method_x_idx[runif(length(method_x_idx)) < base_rate]
  
  # Update their method back to "Others"
  agents@.Data[discontinue_idx, "current_method"] <- "Others"
  
  return(agents)
}

# Run simulation for 5 years
for (year in 2025:2029) {
  agents <- adopt_method_x(agents)
  agents <- discontinue_method_x(agents, high_side_effects = FALSE)
}

# Check final distribution of methods
table(agents@.Data[, "current_method"])
## 
##  Implant Method X     None   Others 
##    29763     4385    14853    50999

This code simulates the adoption and discontinuation of a novel contraceptive technology (Method X) over a 5-year period (2025–2029) using an agent-based model. It models how individuals in a simulated population decide to start using Method X and how some eventually discontinue its use.

Key Actions Taken in the Code:

Define an adoption function (adopt_method_x)

  1. Selects a random percentage (3-7%) of individuals not using any contraceptive method (“None”) and assigns them to Method X each year.

    Define a discontinuation function (discontinue_method_x)

    A default 5% of Method X users stop using it each year.

    If side effects are high, the discontinuation rate doubles to 10%.

  2. Simulate these processes for 5 years (2025–2029)

    Each year, a portion of individuals adopt Method X.

    Some of those already using it decide to discontinue.

  3. Analyze the final contraceptive method distribution

At the end of the 5 years, the code counts how many individuals are still using each method.

**VISUALIZATION OF THE ADOPTION OF METHOD X OVER A 5 YEAR PERIOD**
# Initialize a vector to store the count of Method X users each year
method_x_count <- numeric(length(2025:2029))

# Run simulation for 5 years
for (year in 1:5) {
  agents <- adopt_method_x(agents)
  agents <- discontinue_method_x(agents, high_side_effects = FALSE)
  
  # Count number of agents using "Method X"
  method_x_count[year] <- sum(agents@.Data[, "current_method"] == "Method X")
}

# Plot the adoption trend over 5 years
years <- 2025:2029
plot(years, method_x_count, type="b", col="blue", pch=19, 
     xlab="Year", ylab="Number of Method X Users",
     main="Adoption of Method X Over Time")

# Function to run simulation under different adoption and discontinuation rates
run_sensitivity_analysis <- function(adoption_range, base_discontinuation, high_side_effects) {
  agents_copy <- agents  # Preserve original dataset
  method_x_count <- numeric(length(2025:2029))
  
  for (year in 1:5) {
    # Adjust adoption rate
    adopt_method_x <- function(agents) {
      adoption_rate <- sample(seq(adoption_range[1], adoption_range[2], by=0.01), 1)
      current_methods <- agents@.Data[, "current_method"]
      none_indices <- which(current_methods == "None")
      adopt_indices <- none_indices[runif(length(none_indices)) < adoption_rate]
      agents@.Data[adopt_indices, "current_method"] <- "Method X"
      return(agents)
    }

    # Adjust discontinuation rate
    discontinue_method_x <- function(agents, high_side_effects) {
      discontinuation_rate <- base_discontinuation
      if (high_side_effects) {
        discontinuation_rate <- discontinuation_rate * 2  # Increase if high side effects
      }
      current_methods <- agents@.Data[, "current_method"]
      method_x_indices <- which(current_methods == "Method X")
      discontinue_indices <- method_x_indices[runif(length(method_x_indices)) < discontinuation_rate]
      agents@.Data[discontinue_indices, "current_method"] <- "None"
      return(agents)
    }

    # Run adoption and discontinuation
    agents_copy <- adopt_method_x(agents_copy)
    agents_copy <- discontinue_method_x(agents_copy, high_side_effects)
    
    # Record Method X users
    method_x_count[year] <- sum(agents_copy@.Data[, "current_method"] == "Method X")
  }
  
  return(method_x_count)
}

# Run different scenarios
scenario_1 <- run_sensitivity_analysis(c(0.03, 0.07), 0.05, FALSE)  # Base case
scenario_2 <- run_sensitivity_analysis(c(0.05, 0.10), 0.05, FALSE)  # Higher adoption
scenario_3 <- run_sensitivity_analysis(c(0.03, 0.07), 0.10, FALSE)  # Higher discontinuation
scenario_4 <- run_sensitivity_analysis(c(0.03, 0.07), 0.05, TRUE)   # High side effects

# Plot results
years <- 2025:2029
plot(years, scenario_1, type="b", col="blue", pch=19, ylim=c(0, max(c(scenario_1, scenario_2, scenario_3, scenario_4))), 
     xlab="Year", ylab="Number of Method X Users", main="Sensitivity Analysis: Method X Adoption")
lines(years, scenario_2, type="b", col="green", pch=19)
lines(years, scenario_3, type="b", col="red", pch=19)
lines(years, scenario_4, type="b", col="purple", pch=19)
legend("topright", legend=c("Base Case", "Higher Adoption", "Higher Discontinuation", "High Side Effects"), 
       col=c("blue", "green", "red", "purple"), lty=1, pch=19)