# Chapter 14 # coded by Minato Nakazawa # 3 July 2014. CM2 <- function(a0=18, k=1.0) { ages <- 10:60 em <- fm <- 0.174/k*exp(-4.411*exp(-0.309/k*(ages-a0))) for (i in 2:51) { em[i] <- em[i-1] + (1-em[i-1])*fm[i] } return(list(ages=ages, g=fm, G=em)) } # plot(10:60, CM2(18, 1.0)$G, type="l", ylim=c(0,1), col=1, lty=1, xlab="Age", ylab="Prop. ever married") lines(10:60, CM2(18, 0.5)$G, col=2, lty=2) lines(10:60, CM2(25, 0.5)$G, col=3, lty=3) legend("bottomright", lty=1:3, col=1:3, legend=c("a0=18, k=1.0","a0=18, k=0.5","a0=25, k=0.5")) # 3 parameter CM (p.169 of the text) is included in fmsb package. if (require(fmsb)==FALSE) { install.packages("fmsb"); library(fmsb) } example(CM) # C=1.3 fits 2010 Japan, very strange, because 1.0 means # all women eventually marrying. # example(CT) plot(12:49, ASMFR, type="p") lines(12:49, CT(res[1], res[2]), col="red") # Gompit CPASFR2010 <- cumsum(Jfert$ASFR2010)/sum(Jfert$ASFR2010) GompitCPASFR2010 <- -log(-log(CPASFR2010)) CPASFR1980 <- cumsum(Jfert$ASFR1980)/sum(Jfert$ASFR1980) GompitCPASFR1980 <- -log(-log(CPASFR1980)) layout(t(1:2)) plot(CPASFR1980~CPASFR2010) plot(GompitCPASFR1980~GompitCPASFR2010) # GompitCPASFR2010[35:40] for ages 49-54 are Infinity, so excluded # from regression analysis reg <- lm(GompitCPASFR1980[1:34]~GompitCPASFR2010[1:34]) abline(reg) summary(reg) # predict(reg) = coef(reg)[1] + coef(reg)[2]*GompitCPASFR2010[1:34] (modelCPASFR1980 <- exp(-exp(-predict(reg)))) layout(1) plot(CPASFR1980[1:34] ~ Jfert$Age[1:34], xlab="Age", ylab="Cumulative Proportional ASFRs") lines(Jfert$Age[1:34], modelCPASFR1980, col=2) legend("bottomright", pch=c(1, NA), lty=c(NA, 1), col=1:2, legend=c("Actural 1980 Japan", "Gompit using 2010 std")) Brass <- function(x, c, s) { c*(x-s)*(s+33-x)^2 } Brass(15:49, 1/10000, 12) example("Hadwiger") plot(15:54, Jfert$ASFR2000) lines(15:54, Hadwiger(res[1], res[2], res[3]), col="red")