views of sustainable population mountain

iterative population pyramid.R

library(ggplot2)
library(dplyr)


# start with current year
starting_year <- 2024

# assume this many children born during the lifetime of each woman
children_per_woman <- 2.1

# assume women bear children over this life period
childbearing_ages <- 15:49

# set a random seed
set.seed( starting_year )

# (incorrectly) assume mother's age at childbirth to be uniform
single_year_probability <- children_per_woman / length( childbearing_ages )

# start with a population of one million
total_person_count <- starting_population_size <- 1000000

# (extremely incorectly) start with a uniform age & sex distribution
person_records <-
	data.frame( 
		person_identifier = seq( starting_population_size ) ,
		sex = c( 'M' , 'F' ) ,
		age = sample( 0:119 , starting_population_size , replace = TRUE )
	)

# look at the starting dataset
sapply( person_records , summary )

table( person_records[ , 'sex' ] )

# use life expectancy tables to figure out death probabilities of each individual
# 2021 (from 2024 trustees report)
# https://www.ssa.gov/oact/STATS/table4c6.html
death_probabilities <-
	structure(list(sex = c("M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
	"M", "M", "M", "M", "M", "M", "M", "M", "M", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
	"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F"
	), age = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
	12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
	25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 
	38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 
	51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 
	64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 
	77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 
	90L, 91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 
	102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 
	113L, 114L, 115L, 116L, 117L, 118L, 119L, 120L, 0L, 1L, 2L, 3L, 
	4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 
	18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 
	31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 
	44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 
	57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 
	70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 
	83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L, 95L, 
	96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L, 106L, 
	107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 
	118L, 119L, 120L), prob = c(0.00586, 0.00042, 0.000272, 0.000225, 
	0.000184, 0.000157, 0.00014, 0.000128, 0.000122, 0.000123, 0.000129, 
	0.000138, 0.000164, 0.00022, 0.00031, 0.000446, 0.000637, 0.000868, 
	0.0011, 0.00127, 0.001373, 0.001488, 0.001605, 0.001714, 0.001835, 
	0.001963, 0.002082, 0.002202, 0.00233, 0.002457, 0.002574, 0.002683, 
	0.002787, 0.002881, 0.002974, 0.003074, 0.003175, 0.003295, 0.003444, 
	0.003608, 0.00378, 0.003958, 0.004144, 0.004337, 0.00454, 0.004774, 
	0.005064, 0.005399, 0.005796, 0.006214, 0.006671, 0.007167, 0.007736, 
	0.008351, 0.009035, 0.00977, 0.010567, 0.011398, 0.012291, 0.013224, 
	0.014267, 0.015353, 0.016484, 0.017617, 0.018759, 0.019914, 0.021104, 
	0.022423, 0.023847, 0.025357, 0.02705, 0.02897, 0.031188, 0.033754, 
	0.036747, 0.040563, 0.044308, 0.048498, 0.053229, 0.058778, 0.064617, 
	0.070947, 0.077834, 0.085686, 0.094809, 0.10509, 0.116592, 0.129306, 
	0.142732, 0.157638, 0.174458, 0.193027, 0.21293, 0.232657, 0.251826, 
	0.270943, 0.289756, 0.307998, 0.325393, 0.341662, 0.358746, 0.376683, 
	0.395517, 0.415293, 0.436058, 0.45786, 0.480753, 0.504791, 0.530031, 
	0.556532, 0.584359, 0.613577, 0.644256, 0.676468, 0.710292, 0.745806, 
	0.783097, 0.822251, 0.863364, 0.906532, 1, 0.005063, 0.000393, 
	0.000223, 0.000177, 0.000144, 0.000122, 0.000109, 0.000102, 9.8e-05, 
	9.7e-05, 0.000103, 0.000113, 0.000131, 0.000157, 0.00019, 0.000233, 
	0.000291, 0.000355, 0.000418, 0.000461, 0.000507, 0.000556, 0.00061, 
	0.000666, 0.000722, 0.000775, 0.000831, 0.000889, 0.000952, 0.001025, 
	0.001104, 0.001192, 0.001289, 0.001383, 0.001465, 0.001544, 0.001626, 
	0.001719, 0.001824, 0.00194, 0.002066, 0.002202, 0.002351, 0.002482, 
	0.002622, 0.002789, 0.002994, 0.003219, 0.003467, 0.003729, 0.004011, 
	0.004306, 0.004634, 0.004981, 0.00537, 0.005831, 0.006326, 0.006837, 
	0.007399, 0.008033, 0.008687, 0.009411, 0.010139, 0.010849, 0.01155, 
	0.012216, 0.012952, 0.013844, 0.014863, 0.016028, 0.017329, 0.018859, 
	0.020609, 0.02262, 0.024958, 0.027906, 0.030925, 0.03414, 0.03762, 
	0.041725, 0.046324, 0.051334, 0.056911, 0.063279, 0.070704, 0.079184, 
	0.088697, 0.09924, 0.11048, 0.123078, 0.137152, 0.152605, 0.169494, 
	0.187623, 0.206647, 0.22589, 0.245054, 0.263815, 0.281828, 0.298738, 
	0.316662, 0.335662, 0.355802, 0.37715, 0.399779, 0.423766, 0.449192, 
	0.476143, 0.504712, 0.534994, 0.567094, 0.60112, 0.637187, 0.675418, 
	0.710292, 0.745806, 0.783097, 0.822251, 0.863364, 0.906532, 1
	)), class = "data.frame", row.names = c(NA, -242L))




# start with the current year and go until the year 5,000 if necessary
for( year in seq( starting_year + 1 , 5000 ) ){


	
	childbearing_age_women_count <- 
		nrow( 
			subset( 
				person_records , 
				sex == 'F' & age %in% childbearing_ages
			)
		)

	new_babies_this_year <-
		round( single_year_probability * childbearing_age_women_count , 0 )

	if( new_babies_this_year %% 2 != 0 ) new_babies_this_year <- new_babies_this_year + 1

	# stork delivery
	baby_records <-
		data.frame(
			person_identifier =
				seq( 
					total_person_count + 1 , 
					total_person_count + new_babies_this_year 
				) ,
			sex = c( 'M' , 'F' ) ,
			age = 0
		)
	
	# make everyone one year older
	person_records[ , 'age' ] <- person_records[ , 'age' ] + 1
	
	# combine two datasets
	person_records <- rbind( person_records , baby_records )
	
	# merge on death probabilities
	person_records_with_death_probabilities <-
		merge( person_records , death_probabilities , all.x = TRUE )
		
	# identify decedents
	person_records_with_death_probabilities[ , 'died_this_year' ] <-
		mapply(
			function( a , b ) sample( 0:1 , size = 1 , prob = c( a , b ) ) ,
			1 - person_records_with_death_probabilities[ , 'prob' ] ,
			person_records_with_death_probabilities[ , 'prob' ]
		)
	
	# remove decedents
	person_records <-
		person_records_with_death_probabilities[ 
			person_records_with_death_probabilities[ , 'died_this_year' ] %in% 0 ,
			c( 'age' , 'sex' , 'person_identifier' )
		]

	# aggregate for graphic
	age_sex <- data.frame( table( person_records[ c( 'sex' , 'age' ) ] ) )
	names( age_sex )[ 3 ] <- 'population'
	
	message( paste( year , 'has a population of' , sum( age_sex$population ) ) )

}

# change male population to negative 
age_sex %>%mutate( 
	population = ifelse(sex=="M", population*(-1), 
						population*1))%>% 
	ggplot(aes(x = age,y = population, fill=sex)) +  
	geom_bar(stat = "identity") + 
	coord_flip()+ 
	scale_fill_brewer(type = "seq",palette = 7)+
	theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank()  #remove y axis ticks
        )
	
# readline(prompt="Press [enter] to continue")