# 図5.8 if (!require(fmsb)) { install.packages("fmsb"); library(fmsb) } matplot(Jvital$YEAR, cbind(Jvital$SARPB, Jvital$ACRPB), type="b", pch=c(1, 16), lty=1:2, col=1, xlab="年次", ylab="割合(出産千対)") legend("topleft", pch=c(1, 16), lty=1:2, col=1, legend=c("自然死産率","人工妊娠中絶率")) PrincetonIndex <- function(B, Ws, LB, MWs, ILB, UMWs) { # Hutterite's age class specific fertilities from 15-19 to 45-49. # Source: Coale AJ, Watkins SC: The Decline of Fertility in Europe # Princeton Univ. Press, 1986. pp.154 (Appendix B of Chapter 2) HFs <- c(0.300, 0.550, 0.502, 0.447, 0.406, 0.222, 0.061) If <- B/sum(Ws*HFs) Ig <- LB/sum(MWs*HFs) Ih <- ILB/sum(UMWs*HFs) Im <- sum(MWs*HFs)/sum(Ws*HFs) return(list(If=If, Ig=Ig, Ih=Ih, Im=Im)) } # Example for Princeton Index # Wachter KW "Essential Demographic Methods" # Chapter 6, Section 7 "Princeton Indices" # Hutterites' fertility for ages 30-34 is wrongly given as 0.407 # The data shows Berlin in 1900 B.Ws <- c(91358, 114464, 99644, 88886, 75729, 66448, 54485) B.MWs <- c(1538, 28710, 55417, 62076, 55293, 47197, 36906) B.B <- 49638 B.LB <- 42186 PrincetonIndex(B.B, B.Ws, B.LB, B.MWs, B.B-B.LB, B.Ws-B.MWs) # Coale-Trussellモデルの当てはめ Y <- sprintf("ASMFR%d",1950+0:13*5) Mm <- data.frame( M=numeric(14), m=numeric(14)) for (i in 1:14) { ASMFR <- c(0, 0, 0, Jfert[1:35, Y[i]]) res <- fitCT(, ASMFR) FLAG <- res[4] while (FLAG>0) { res <- fitCT(res[1:2], ASMFR) FLAG <- res[4] } Mm$M[i] <- res[1] Mm$m[i] <- res[2] } rownames(Mm) <- Y print(Mm) # HadwigerモデルのASMFRへの当てはめ Year <- 1950+0:13*5 Years <- sprintf("ASMFR%d", Year) Hadpars <- data.frame(Year = Year, a = numeric(length(Year)), b = numeric(length(Year)), c = numeric(length(Year)), RMSE = numeric(length(Year))) for (i in 1:length(Year)) { res <- fitHad(,Jfert[,Years[i]]) while (res[5]>0) { res <- fitHad(,Jfert[,Years[i]]) } Hadpars[i,] <- c(Year[i], res[1], res[2], res[3], res[4]) } spr <- function(X) { sprintf("%d & %5.3f & %5.3f & %5.3f & %5.4f\\\\\n", X[1], X[2], X[3], X[4], X[5])} cat(apply(Hadpars, 1, spr)) # Hadwigerモデルの日本のASFRへの当てはめ Years2 <- sprintf("ASFR%d", Year) Hadpars2 <- data.frame(Year = Year, a = numeric(length(Year)), b = numeric(length(Year)), c = numeric(length(Year)), RMSE = numeric(length(Year))) for (i in 1:length(Year)) { res <- fitHad(, Jfert[,Years2[i]]) while (res[5]>0) { res <- fitHad(, Jfert[,Years2[i]]) } Hadpars2[i,] <- c(Year[i], res[1], res[2], res[3], res[4]) } cat(apply(Hadpars2, 1, spr))