Demography Package In R Download
Using the demography package
This package from Rob J Hyndman from Monash University in Australia includes functions for common calculations used in demography such as lifetable calculations; Lee-Carter modelling and variants; functional data analysis of mortality rates, fertility rates, net migration numbers; and stochastic population forecasting. The documentation for the package discusses the functions in alphabetical order using the sample datasets provided. Functions for preparing your data are not at the beginning of the alphabet, so it's easy to overlook them and not know how to make use of the package with real data. This document will help you isolate the functions that load the data and talk a bit about the data structure required to use the package.
Play with the sample data provided in the demography package
require(demography) plot(fr.mort)
plot(fr.mort, series = "male")
The demogdata object
The sample data included in the demography package (aus.fert, fr.mort, fr.sm) are of the class "demogdata". This is a complex, custom data structure that facilitates the use of the functions in the package. If we take a look at the structure fr.mort data we see:
str(fr.mort) ## List of 7 ## $ type : chr "mortality" ## $ label : chr "FRATNP" ## $ lambda: num 0 ## $ year : int [1:191] 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 ... ## $ age : num [1:111] 0 1 2 3 4 5 6 7 8 9 ... ## $ rate :List of 3 ## ..$ total : num [1:111, 1:191] 0.2053 0.0467 0.0341 0.023 0.016 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## ..$ female: num [1:111, 1:191] 0.187 0.0467 0.0339 0.0229 0.016 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## ..$ male : num [1:111, 1:191] 0.2229 0.0467 0.0343 0.0232 0.0161 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## $ pop :List of 3 ## ..$ total : num [1:111, 1:191] 834355 782273 714855 686823 674203 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## ..$ female: num [1:111, 1:191] 408224 382452 351454 337733 331576 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## ..$ male : num [1:111, 1:191] 426130 399821 363401 349090 342627 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:111] "0" "1" "2" "3" ... ## .. .. ..$ : chr [1:191] "1816" "1817" "1818" "1819" ... ## - attr(*, "class")= chr "demogdata" This is complicated! But you don't have to build all these objects, you just need to understand the simple, rectangular structure of the data files in the Human Mortality database (see below) then use one of the two functions in the package to convert the data into the structure above.
Specifying different plot types for demogdata
par(mfrow = c(2, 2)) plot(aus.fert, plot.type = "time") plot(aus.fert, plot.type = "functions") plot(sex.ratio(fr.mort), ylab = "Sex ratios (M/F)")
Plotting subsets
par(mfrow = c(1, 2)) plot(fr.mort, year = 1917:1945, series = "female") plot(fr.mort, year = 1917:1945, series = "male")
Misc snippets
# LIFETABLES france.lt <- lifetable(fr.mort) plot(france.lt)
# Print some of the life table lt1990 <- print(lifetable(fr.mort, year = 1990)) ## Period lifetable for FRATNP : total ## ## Year: 1990 ## mx qx lx dx Lx Tx ex ## 0 0.0075 0.0074 1.0000 0.0074 0.9931 76.8507 76.851 ## 1 0.0006 0.0006 0.9926 0.0006 0.9923 75.8576 76.424 ## 2 0.0004 0.0004 0.9920 0.0004 0.9918 74.8653 75.469 ## 3 0.0003 0.0003 0.9916 0.0003 0.9915 73.8735 74.500 ## 4 0.0002 0.0002 0.9913 0.0002 0.9912 72.8821 73.521 ## 5 0.0002 0.0002 0.9911 0.0002 0.9910 71.8909 72.538 ## 6 0.0002 0.0002 0.9909 0.0002 0.9908 70.8999 71.553 ## 7 0.0002 0.0002 0.9907 0.0002 0.9906 69.9091 70.568 ## 8 0.0002 0.0002 0.9905 0.0002 0.9904 68.9185 69.580 ## 9 0.0002 0.0002 0.9903 0.0002 0.9902 67.9281 68.591 ## 10 0.0002 0.0002 0.9902 0.0002 0.9901 66.9379 67.604 ## 11 0.0002 0.0002 0.9900 0.0002 0.9899 65.9478 66.617 ## 12 0.0002 0.0002 0.9898 0.0002 0.9897 64.9580 65.628 ## 13 0.0002 0.0002 0.9896 0.0002 0.9895 63.9683 64.640 ## 14 0.0002 0.0002 0.9894 0.0002 0.9893 62.9787 63.652 ## 15 0.0004 0.0004 0.9892 0.0004 0.9890 61.9895 62.668 ## 16 0.0004 0.0004 0.9888 0.0004 0.9886 61.0005 61.691 ## 17 0.0006 0.0006 0.9884 0.0006 0.9881 60.0119 60.718 ## 18 0.0008 0.0008 0.9878 0.0008 0.9874 59.0238 59.755 ## 19 0.0009 0.0009 0.9870 0.0009 0.9866 58.0364 58.801 ## 20 0.0009 0.0009 0.9861 0.0009 0.9857 57.0498 57.852 ## 21 0.0010 0.0010 0.9853 0.0010 0.9848 56.0642 56.903 ## 22 0.0010 0.0010 0.9843 0.0010 0.9838 55.0794 55.960 ## 23 0.0010 0.0010 0.9833 0.0010 0.9828 54.0957 55.017 ## 24 0.0010 0.0010 0.9823 0.0010 0.9818 53.1129 54.073 ## 25 0.0010 0.0010 0.9813 0.0010 0.9807 52.1311 53.127 ## 26 0.0011 0.0011 0.9802 0.0011 0.9797 51.1504 52.181 ## 27 0.0011 0.0011 0.9792 0.0011 0.9786 50.1707 51.237 ## 28 0.0011 0.0011 0.9781 0.0010 0.9776 49.1920 50.294 ## 29 0.0012 0.0012 0.9771 0.0011 0.9765 48.2145 49.347 ## 30 0.0012 0.0012 0.9759 0.0011 0.9754 47.2380 48.404 ## 31 0.0013 0.0013 0.9748 0.0012 0.9742 46.2626 47.459 ## 32 0.0012 0.0012 0.9735 0.0012 0.9730 45.2885 46.519 ## 33 0.0014 0.0014 0.9724 0.0013 0.9717 44.3155 45.575 ## 34 0.0014 0.0014 0.9710 0.0014 0.9703 43.3438 44.636 ## 35 0.0015 0.0015 0.9696 0.0014 0.9689 42.3735 43.700 ## 36 0.0016 0.0016 0.9682 0.0016 0.9674 41.4045 42.764 ## 37 0.0017 0.0017 0.9666 0.0016 0.9658 40.4371 41.833 ## 38 0.0019 0.0019 0.9650 0.0019 0.9641 39.4713 40.902 ## 39 0.0019 0.0019 0.9632 0.0018 0.9623 38.5072 39.980 ## 40 0.0021 0.0021 0.9614 0.0020 0.9604 37.5449 39.054 ## 41 0.0022 0.0022 0.9594 0.0021 0.9583 36.5846 38.134 ## 42 0.0024 0.0024 0.9572 0.0023 0.9561 35.6263 37.218 ## 43 0.0027 0.0026 0.9549 0.0025 0.9537 34.6702 36.306 ## 44 0.0029 0.0029 0.9524 0.0028 0.9510 33.7165 35.401 ## 45 0.0031 0.0031 0.9496 0.0030 0.9482 32.7655 34.503 ## 46 0.0034 0.0034 0.9467 0.0033 0.9450 31.8173 33.609 ## 47 0.0036 0.0036 0.9434 0.0034 0.9417 30.8723 32.724 ## 48 0.0039 0.0039 0.9400 0.0036 0.9382 29.9305 31.840 ## 49 0.0038 0.0038 0.9364 0.0036 0.9346 28.9923 30.962 ## 50 0.0048 0.0048 0.9328 0.0044 0.9306 28.0578 30.080 ## 51 0.0049 0.0049 0.9283 0.0045 0.9261 27.1272 29.221 ## 52 0.0055 0.0054 0.9238 0.0050 0.9213 26.2011 28.362 ## 53 0.0059 0.0059 0.9188 0.0054 0.9161 25.2798 27.514 ## 54 0.0064 0.0063 0.9134 0.0058 0.9105 24.3637 26.675 ## 55 0.0072 0.0072 0.9076 0.0065 0.9043 23.4533 25.842 ## 56 0.0074 0.0074 0.9011 0.0067 0.8977 22.5490 25.025 ## 57 0.0079 0.0079 0.8944 0.0070 0.8909 21.6513 24.208 ## 58 0.0089 0.0088 0.8873 0.0078 0.8834 20.7604 23.396 ## 59 0.0096 0.0095 0.8795 0.0084 0.8753 19.8770 22.600 ## 60 0.0101 0.0101 0.8711 0.0088 0.8667 19.0016 21.812 ## 61 0.0110 0.0109 0.8624 0.0094 0.8577 18.1349 21.029 ## 62 0.0119 0.0119 0.8529 0.0101 0.8479 17.2772 20.256 ## 63 0.0128 0.0127 0.8428 0.0107 0.8375 16.4294 19.493 ## 64 0.0136 0.0135 0.8321 0.0112 0.8265 15.5919 18.738 ## 65 0.0143 0.0142 0.8209 0.0116 0.8151 14.7654 17.987 ## 66 0.0154 0.0152 0.8092 0.0123 0.8031 13.9503 17.239 ## 67 0.0166 0.0165 0.7969 0.0131 0.7903 13.1472 16.498 ## 68 0.0182 0.0181 0.7838 0.0142 0.7767 12.3569 15.766 ## 69 0.0194 0.0192 0.7696 0.0148 0.7622 11.5802 15.047 ## 70 0.0227 0.0224 0.7548 0.0169 0.7463 10.8180 14.332 ## 71 0.0225 0.0222 0.7379 0.0164 0.7297 10.0717 13.649 ## 72 0.0268 0.0265 0.7215 0.0191 0.7120 9.3420 12.948 ## 73 0.0290 0.0285 0.7024 0.0200 0.6924 8.6300 12.286 ## 74 0.0292 0.0288 0.6824 0.0196 0.6725 7.9376 11.633 ## 75 0.0377 0.0370 0.6627 0.0245 0.6505 7.2651 10.963 ## 76 0.0387 0.0379 0.6382 0.0242 0.6261 6.6146 10.364 ## 77 0.0440 0.0430 0.6140 0.0264 0.6008 5.9885 9.754 ## 78 0.0481 0.0470 0.5876 0.0276 0.5738 5.3878 9.170 ## 79 0.0545 0.0530 0.5599 0.0297 0.5451 4.8140 8.597 ## 80 0.0610 0.0592 0.5302 0.0314 0.5145 4.2689 8.051 ## 81 0.0695 0.0672 0.4989 0.0335 0.4821 3.7544 7.526 ## 82 0.0785 0.0756 0.4653 0.0352 0.4478 3.2723 7.032 ## 83 0.0880 0.0843 0.4302 0.0362 0.4121 2.8245 6.566 ## 84 0.0986 0.0939 0.3939 0.0370 0.3754 2.4125 6.124 ## 85 0.1093 0.1036 0.3569 0.0370 0.3384 2.0370 5.707 ## 86 0.1234 0.1163 0.3200 0.0372 0.3014 1.6986 5.309 ## 87 0.1369 0.1281 0.2828 0.0362 0.2646 1.3972 4.942 ## 88 0.1515 0.1408 0.2465 0.0347 0.2292 1.1326 4.594 ## 89 0.1717 0.1582 0.2118 0.0335 0.1951 0.9034 4.265 ## 90 0.1890 0.1727 0.1783 0.0308 0.1629 0.7084 3.973 ## 91 0.2068 0.1874 0.1475 0.0276 0.1337 0.5455 3.698 ## 92 0.2336 0.2092 0.1199 0.0251 0.1073 0.4118 3.435 ## 93 0.2563 0.2272 0.0948 0.0215 0.0840 0.3044 3.211 ## 94 0.2833 0.2482 0.0733 0.0182 0.0642 0.2204 3.008 ## 95 0.3009 0.2615 0.0551 0.0144 0.0479 0.1562 2.836 ## 96 0.3257 0.2801 0.0407 0.0114 0.0350 0.1083 2.664 ## 97 0.3529 0.3000 0.0293 0.0088 0.0249 0.0734 2.506 ## 98 0.3618 0.3064 0.0205 0.0063 0.0174 0.0485 2.365 ## 99 0.4087 0.3394 0.0142 0.0048 0.0118 0.0311 2.189 ## 100 0.4862 1.0000 0.0094 0.0094 0.0193 0.0193 2.057 # Lee-Carter model france.LC <- lca(fr.mort) france.fcast <- forecast(france.LC) france.lt.f <- lifetable(france.fcast) plot(france.lt.f) 
Forecasting
### Functional Model Forecast Birth cohort lifetables, 1900-1910 france.clt <- lifetable(fr.mort, type = "cohort", age = 0, years = 1900:1910) # Model Forecasting france.fit <- fdm(fr.mort, order = 2) france.fcast <- forecast(france.fit, 50) plot(france.fcast) models(france.fcast) ## ## -- Coefficient 1 -- ## Series: xx[, i] ## ARIMA(1,1,1) with drift ## ## Coefficients: ## ar1 ma1 drift ## 0.616 -0.782 -0.127 ## s.e. 0.178 0.141 0.043 ## ## sigma^2 estimated as 1.09: log likelihood=-275.9 ## AIC=559.9 AICc=560.1 BIC=572.9 ## ## -- Coefficient 2 -- ## Series: xx[, i] ## ARIMA(1,1,1) ## ## Coefficients: ## ar1 ma1 ## 0.741 -0.951 ## s.e. 0.071 0.035 ## ## sigma^2 estimated as 0.448: log likelihood=-193.6 ## AIC=393.3 AICc=393.4 BIC=403 # Simulation france.sim <- simulate(france.fcast, nsim = 100) # Plot.fmforecast plot(france.fcast) 
plot(france.fcast, "c") 
plot(france.fcast, "v") 
Using the Human Mortality Database
The hmd.mx function makes it easy to read data directly from the Human Mortality Database website http://www.mortality.org/ without even opening your browser and turn them into demogdata objects. The function can read the appropriate single-year interval period files ("1x1"s) for mortality, fertility and migration, you just need to specify the country code.
You will need to register for access to the Human Mortality Database. It's free. Then you can use R code to pull data and analyze. Use this code below to try it.
# Read more data from Human Mortality Database (substitute your HMD id/pw for x's) poland<-hmd.mx("POL", "yourHMDemail@psu.edu", "yourHMDpassword", "Poland") Reading data from a text file
The HMD has a limited number of countries. If you have figures for death rates, exposure to risk, population, etc from other countries not in the database, just record them in a tab-delimited file. You can load the numbers from that file into a demogdata object using the read.demogdata function. You need to setup your text files in a particular format for the function and that format is based on the HMD raw text files. Before you prepare your own data, it is helpful to download a couple of text files from the HMD and take a look at them. So let's make another copy of the Polish data but not with the hmd.mx function. Instead, open a browser and go to the HMD Poland page and download two Period files for single age and single year intervals(1x1's): Exposure to Risk 1x1(Exposures_1x1.txt) and Death Rates 1x1(Mx_1x1.txt) tables.
# To demonstrate how to read data from flat files, download the raw tab-delimited files from HMD poland2<-read.demogdata("Mx_1x1.txt", "Exposures_1x1.txt", type="mortality", label="Poland2") The first few lines of the Death rates file look like this:
HMD rectangular file structure for death rate data Poland, Death rates (period 1x1) Last modified: 12-Nov-2010, MPv5 (May07) Year Age Female Male Total 1958 0 0.065861 0.083974 0.075143 1958 1 0.004744 0.005040 0.004895 1958 2 0.001737 0.002037 0.001890 1958 3 0.001192 0.001463 0.001330 1958 4 0.000853 0.001073 0.000965 1958 5 0.000708 0.000943 0.000828 1958 6 0.000599 0.000850 0.000727 # Compare par(mfrow=c(2,1)) plot(poland) plot(poland2) If you plot(poland) and plot(poland2) you will see the same results.
Now you can use the Mx_1x1.txt and Exposures_1x1.txt files as templates to enter your own data with a text editor, Excel or to create an extraction from a SAS or Stata dataset. For more information on how to setup your text file for this function, refer to the HMD Explanatory Notes and the R help on the ?read.demogdata function.
Want to see what code is behind these functions? Type functions w/o parms.
lifetable plot.demogdata fdm Source: http://rstudio-pubs-static.s3.amazonaws.com/10379_1a1858a1259b4c019fb1092922e307b3.html
Posted by: naturalmils.blogspot.com