mpascariu / mortalitylaws Goto Github PK
View Code? Open in Web Editor NEWFit and compare the most popular human mortality laws - R package
License: Other
Fit and compare the most popular human mortality laws - R package
License: Other
On 14 Dec 2018, at 14:57, Patrick G. wrote:
Dear Marius - I have been testing MortalityLaws for some extensions on single age life tables, and I noted that your LifeTable function uses a strong assumption on ax for age 0. Right now your assumption using a scalar works fine for age 1 onward, but for age 0 you should allow the user either to input their a0 value or use the CD-West assumption or Andreev variation.
Hi Marius,
I've noticed that the package will use Heligman-Pollard to generate the 8 parameters and then generate the ASDR for every ages by substituting back the parameter values with x = age level.
However, I'm currently modeling ASDR according to age-groups (0, 1-4, 5-9, 10-14, 15-19 etc). And thus, when using the package I would have to rename the age column to the first age of every age group (0, 1, 5, 10 etc) so that R would be able to read. I'm wondering I should make some modification to the codes if I'm modeling ASDR according to age-groups or just using the package as it is would be fine, since it might seem that the fitted values generated would be the mortality rate for the age levels instead of the age groups.
Also, I tried to generate the AIC and BIC criterion with my current work but it returns 'NA'. I'm not sure what's the issue
Here's the code, might to be a long-winded procedure, it's because I'm trying to be as close to how the data are framed in the Introduction manual provided:
Appreciate your input on the matter. Thanks!
Hi!
I got a question about gompertz's formula:
You wrote:
"
Table 1: Parametric functions build in the MortalityLaws package
Mortality lawsPredictor
Gompertz
µ(x) = AeBx
Gompertz
Code
µ(x) = σ1 exp
"
Did you mean that gompertz formula is: µ(x) = A * B^x (Parameter "A" multiplied parameter "B" raised to "X" (X being age) , or would you define it another way ?
I also have doubts about how the formulas of Makeham and Thiele are for you, and based on what formulas are made the codes that return the parameters.
BTW, thanks for building and share a package about biometrics!
I hope you can read me soon,
Sheila
When using MortalityLaw() and then fitted(), the values are not corresponding to the true value. This also affect LifeTable()
library(MortalityLaws)
inf = 7
sup = 9
x = c(5:10)
qx = c(0.000436, 0.000438, 0.00044, 0.000443, 0.000446, 0.000449)
ml1 <- MortalityLaw(x, qx = qx, law = 'gompertz', fit.this.x = c(inf:sup))
ml1
coef(ml1)
fitted(ml1)
In the code I put values corresponding to ages 5 to 10, and fitting only 7 to 9 with Gompertz Law.
But when calling fitted() the corresponding values are not correct, the display values are 0.0004341233 for 5 (this must be 0 instead of 5), 0.0004370643 for 6 (this must be 1 instead of 6), ...
It seems that instead of using x=5 (and the respective value), it starts in x=0. This could not be a problem, but it also affect LifeTable() making the ages to be displaced.
Using the coefficients and the formula for Gompertz, with x=5 the corresponding mu[x]=Aexp[Bx] value must be 0.00045207065
hello,dr.mpascariu, i am a learner from peking university. i tried to learn how to use R to correct the mortality.i am very glad to find you page,but a little confused to use the code of motalitylaws,because there are so many files,how to use it? i am looking forward for your reply,thank you
The actual version of the package uses a simple method to close a life table. It works well in the most common cases but when dealing with abridged life-tables and various indices more particular methods and flexibility is required.
In order to come up with a solution, we have to figure out:
This issue is related to issue #7 .
Hi, I'm facing some issues when fitting and forecasting using the HP model as I'm still new to using R:
For reference, I'm looking for something like this when plotting the HP parameters:
Appreciate your help on the matter, thanks in advance :)
There are some places where we determine action based on class:
if (K$iclass == "numeric") {
inside MortalityLaw()
, where the value of iclass
comes from find.my.case()
. The below code comes from a DemoTools
example that previously didn't throw warnings. Now class of mx_matrix
will return as c("matrix" ,"array" )
, so we have a length warning in the above check. However, the following situations might also not behave as desired. 1) is a normal case. 2) used to not throw warnings, but now does. 3) throws no warning, but might not be handled properly. 4) warns and might not be handled properly. I'm not sure about the potential mishandling, though, as perhaps those calls are picked up elsewhere or downstream. If you tell me a potential better strategy I may be able to make a usable PR. So far my PR makes case (3) succeed where it used to error, but these appear to be related loose ends.
mx1 <- c(.0070, .0082, .0091, .0096, .0108, .0122, .0141, .0150, .0165, .0186,
.0205, .0229, .0259, .0294, .0334, .0379, .0426, .0482, .0550, .0628,
.0716, .0806, .0897, .1003, .1149, .1264, .1558, .1563, .1812, .2084,
.2298, .2536, .2813, .3143, .3352, .3651, .4128)
x1 <- 65:101
# 1 works fine
MortalityLaw(x =x1, mx = mx1, law = 'makeham')
# 2) try for matrix of inputs
mx_matrix <- matrix(rep(mx1, 3), ncol = 3) %*% diag(c(1, 1.05, 1.1))
dimnames(mx_matrix) <- list(age = x1, year = c("year1", "year2", "year3"))
MortalityLaw(x =x1, mx = mx_matrix, law = 'makeham') # warns
# 3)
# no warning, BUT, we have something that should be treated as a numeric vector, but class evaluates as array, so it won't execute that code. dput() will create output like this sometimes...or always?
mx1.1 <- structure(mx1, .Dim = length(mx1))
MortalityLaw(x =x1, mx = mx1.1 , law = 'makeham')
# 4)
# 1-column matrix should be treated as numeric vector, no?, but here we get the same warning as the 2d case above
MortalityLaw(x =x1, mx =as.matrix(mx1) , law = 'makeham')
Below is the code snippet in questions from LifeTable.R of the LifeTable.core method:
Tx <- rev(cumsum(rev(Lx)))
ex <- Tx/(lx - dx*(ax/nx))
Here Tx is the sum of person-years lived from the current row on down and ex should be the expected total life years lived for a person starting that row. It looks like you are dividing by the current row's person-years lived instead of the number of people who started the interval.
US life tables use the people starting the interval as the denominator. I apologize if I'm missing something.
Will you be updating/resubmitting to CRAN in future?
Solution needed for checking the available countries/populations in HMD.
If a new country is added in the database availableHMD() function fails to find it, because I haven't found an elegant method to read all the country codes.
No problem with the codes that I already know that are there.
Let me know if you have a solution. Thanks.
Marius
Hello, I'm trying to obtain the fitted values based on the parameters that the heligman-Pollard model give me but i can't obtain the qx that R printed. Those are my parameters and the values that I calculated and the ones that the program printed.
parameter Age calculated values Model Values
A 0.000056232 0 #¡NUM! 0.000124871
B 0.006158410 1 0.000056224 0.000094944
C 0.009473579 2 0.000052732 0.000094705
D 0.000624343 3 0.000050779 0.000096278
E 16.407580000 4 0.000049450 0.000098754
F 25.115760000 5 0.000048607 0.000101875
G 0.000035750
H 1.084004000
I hope you can help me with this issue
Consider including the Strehler-Mildvan model (1960) in the package. The model is discussed in Finkelstein (2012).
The conversion from mx to qx is done using the assumption that mx is equal to a constant force of mortality in the age interval. This gives:
nqx = 1 - exp(-n*nMx)
This implies that the separation factor, nax, is:
nax = n + 1/nMx - n/nqx
However, the 1a0 and 4a1 are replaced by the Coale-Demeny (see p. 20 of Coale and Demeny, 1983) west formulas (here for males):
1a0 = 0.33 if q0 >= 0.100
= 0.045 + 2.875 * q0
This seems to have been transformed to make it a function of m0 as:
1a0 = 0.33 if 1m0 >= 0.107
= 0.045 + 2.684 * 1m0
This does not give the correct results. I am working on the revised formula. Should be ready tomorrow.
Dear Mr. Pascariu:
Despite the registration in the new HMD site, when using the package MortalityLaws is not possible to import any dataset from any country. We following error appears:
Ex data for FRATNP state in the 1x1 format was not to be found. We have been looking here:
https://former.mortality.org/hmd/ESP/Exposures_1x1.txt
I am trying to simulate for other countries but I cannot import the datasets.
In addition, I was checking the MortalityEstimate package and would like to adapt it to different countries instead of England and Wales, France, Sweden and USA. Based on that package I am trying to modify the functions just to import the country data I am trying to analyze but I got the same error when using the MortalityLaws package.
Any ideas on what could be the problem?
Thanks in advance for your time and support.
Hi,
I'm working on Makeham's law and I'm having trouble getting the standard errors of the parameters, being the "stats" component the one that returns the information of interest. This is the only component of the list that cannot be obtained as part of the output. How can I get the standard error?
Here's the code I'm working with:
'''
mort <- structure(list(Age = 27:50, qx_obs = c(0.000224299, 0.000259561,
0, 0, 6.59e-05, 0.000166842, 0.000165868, 0.00010441,
0.000183992, 0.000165593, 0.00012395, 8.64e-05, 0.000206636,
0.000294348, 0.000288283, 0.000287975, 0.000278501,
0.000341287, 0.000422662, 0.00049801, 0.000478205,
0.00067597, 0.000590254, 0.00070058)),
row.names = c(NA, 24L), class = "data.frame")
Makeham<- MortalityLaws::MortalityLaw(x = mort$Age, qx = mort$qx_obs, law = 'makeham', opt.method = 'LF2')
summary(Makeham)
ls(Makeham) #The "stats" component is not part of the list.
'''
Should the a0 assumption not default to "total" in the case that no sex is specified? As is, the result is rather close to a midpoint in the interval.
Here are some example data:
library(MortalityLaws)
mx <- c(0.00439, 0.00016, 0.00009, 0.00012, 0.00023, 0.00031, 0.00036,
0.00047, 0.00065, 0.00095, 0.00148, 0.00243, 0.00374, 0.0058,
0.00893, 0.01444, 0.02383, 0.04183, 0.07672, 0.14308, 0.2426,
0.37998, 0.53975, 0.67115)
x <- c(0,1,seq(5,110,5))
# no sex specified:
LT1 <- LifeTable(mx=LT$mx,x=LT$Age)
# sex given as total:
LT2 <- LifeTable(mx=LT$mx,x=LT$Age, sex = "total")
LT1$lt$ax[1]
LT2$lt$ax[1]
Thanks! Alyson & Tim
The convertFx
work only when a table is provided in input (2 or more years). It should be able to return results for vector inputs as well.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.