Hi Prof. Hyndman

I want to raise a concern regarding the fitted values (and subsequent residuals) calculated by demography::lca. If chooseperiod is set to TRUE in the original function call then the fitted values of the output will always be in terms of log death rates - regardless of the restype (response type) parameter specified.

A snippet from the relevant R code below illustrates the issue. If chooseperiod is set to TRUE then the function uses the specified break method to identify the appropriate fitting period - it returns a LCA model of the corresponding period but with chooseperiod set to FALSE. Under both break methods, the response type parameter is not inherited from the original function call and will therefore assume its default value - log death rates. This creates an inconsistency in the nature of the fitted values (and residuals) returned versus those specified in the model. This is particularly problematic if users of this package want to compare in-sample values between various models - not in terms of log death rates.

   if(chooseperiod)
    {
        if(breakmethod=="bai")
        {
            x <- 1:m
            # Find breakpoints
            bp <- strucchange::breakpoints(kt ~ x)$breakpoints
            # Omit breakpoints less than minperiod from end
            bp <- bp[bp <= (m-minperiod)]
            bestbreak <- max(bp)
            return(lca(data,series,year[(bestbreak+1):m],ages=ages,max.age=max.age,
                adjust=adjust,interpolate=interpolate,chooseperiod=FALSE,scale=scale))
        }
        else
        {
            RS <- devlin <- devadd <- numeric(m-2)
            for(i in 1:(m-2))
            {
                tmp <- lca(data,series,year[i:m],ages=ages,max.age=max.age,adjust=adjust,chooseperiod=FALSE,interpolate=interpolate,scale=scale)
                devlin[i] <- tmp$mdev[2]
                devadd[i] <- tmp$mdev[1]
                RS[i] <- (tmp$mdev[2]/tmp$mdev[1])
            }
            bestbreak <- order(RS[1:(m-minperiod)])[1]-1
            out <- lca(data,series,year[(bestbreak+1):m],ages=ages,max.age=max.age,
                adjust=adjust,chooseperiod=FALSE,interpolate=interpolate,scale=scale)
            out$mdevs <- ts(cbind(devlin,devadd,RS),start=startyear,deltat=deltat)
            dimnames(out$mdevs)[[2]] <- c("Mean deviance total","Mean deviance base","Mean deviance ratio")
            return(out)
        }
    }

Fortunately, this issue only affects in-sample values. All of the forecasting and comparison functions within the package will be based on out-sample data and ignore the in-sample residuals. Furthermore, the solution would appear to be as simple as adding the response type parameter to the models returned - so that it is inherited from the original function call.

This is the first instance of this issue that I encountered and hopefully the last. However, I would recommend revising the parameter specification of all recursive functions so that they are consistent with the original function call.

I hope I explained my concern well enough. Please let me know if I can provide any additional information to clarify my concern.

Regards

Can I apply lca to abridged mortality rate? Age: 0, 1, 5, ...100 Years: 5 year steps. I tried it with 10 period data and following is what I get when I run the following codes

data <- read.demogdata(file="rate.txt", popfile="popfile.txt", type="mortality", label="NepalMale"
, max.mx = 10, skip = 0, popskip = 0, lambda=0, scale=1000)

nepalmale <- lca(data=data, series=names(data$rate)[1], years=data$year,
ages=data$age, max.age=100,
adjust = c("dt", "dxt", "e0", "none")[1], chooseperiod=FALSE,
minperiod=9, breakmethod=c("bai","bms"), scale = FALSE, restype = c("logrates", "rates", "deaths")[2],
interpolate = FALSE)

Error in matrix(mx, nrow = n, ncol = m) :
invalid 'ncol' value (too large or NA)

Hi @robjhyndman
I'm currently pushing HMDHFDplus version 2.0.3 to CRAN, it should be up in a day or two or three. I suggest you change the DESCRIPTION file text HMDHFDplus (>= 2.0.0) to HMDHFDplus (>= 2.0.3) once it's up on CRAN. The main fix is to include national subpopulations that weren't getting scraped in the previous update. e.g. Scotland.

I am trying to implement the lee-carter model for fertility but the lca code is limited only to mortality data.

I tried using fdm(order=1, lambda=Inf, method="M") as suggested by Hyndman and Ullah (2007, p.4953).

I tested this with the france mortality data provided in the package but this does not provide same result as lca.

Hello, i have some models from the StMoMo package and i need to make lifetables from them. When i use the function lifetable i get the error that the data are not demogdata. Can you please help me?

My code is : ```
ages.fit<-55:89
wxt<- genWeightMat(ages = 55:89, years = grstomo.male$years, clip = 3)
rh<-rh(link = "log", cohortAgeFun = "1")
rhfit.male<- fit(rh, data = grstomo.male, ages.fit = ages.fit, wxt=wxt)

When using the code smooth.demogdata and setting weight="TRUE", I am having the following error:

total.sm=smooth.demogdata(total, method="mspline",age.grid=0:85,weight="TRUE",b=60)

Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) :
A term has fewer unique covariate combinations than specified maximum degrees of freedom

How to interpret this error and how to solve it?

Hello,
I encountered an error when I tried to forecast an fdm object with Japan mortality rates.

The code:

jpn <- hmd.mx("JPN", "username", "password", "Japan")
jpn2 <- extract.years(jpn, 1947:2004)
jpn2 <- extract.ages(jpn2, 0:100, combine.upper = TRUE)
fit <- fdm(jpn2, series = "total", method = "M")
fcast <- forecast(fit, 11)

which would return
Error in fitted[, i] <- fitted(barima) :
number of items to replace is not a multiple of replacement length

The same error was also encountered when using mortality rates of Austria, Belgium, Bulgaria, Hungary, Ireland and Spain that were cleaned or transformed to avoid 0 and NA values.
However, I did not receive an error when using mortality rates of Australia, Italy, France and Scotland.

Any help to address these inconsistencies would be much appreciated. Thank you in advance.

Hi Prof Hyndman, I noticed that NAs are introduced when read in the dataset using hmd method, it might because the HMD dataset now includes an age group which is named "110+", and when converting the 'age' column into a numeric vector, "110+" is converted into NA by coercion.

Hi Rob,

I'm facing another issue where I can't seem to generate the AIC & BIC criterion due to object class.
Here's the error message:

Should I calculate the criterion using another package? As it seems like other object classes (StMoMo, MortalityLaw) are able to generate AIC and BIC values.

Appreciate your input on the matter. Thanks!

Hi Rob,

I am using the fdm model to forecast age specific death rates. I have noticed an interesting behaviour of the forecast trajectory in the first 3 years, in the sense that the slope changes significantly from year to year. At some ages the change is more pronounced. See an example below for age 0 and 95.

Could you please indicate if this is normal and what is the reason behind it?
The jumpchoice doesn't seem to influence it. Also this behaviour is not seen in the lca model.

Thanks.
Marius

The figures can be reproduced using this code:

# Thu Aug 30 11:10:36 2018 --------- Marius D. Pascariu ---
remove(list = ls())
library(demography)

M <- fdm(fr.mort)
P <- forecast(france.fit, h = 5, jumpchoice = "actual")

x = 0
plot(P$rate$total[x + 1, ], type = "l", 
     ylab = paste0("m[", x,"]"), 
     xlab = "Forecast Year",
     main = paste("Forecast: Age-Specific Death Rate at age", x))

x = 95
plot(P$rate$total[x + 1, ], type = "l", 
     ylab = paste0("m[", x,"]"), 
     xlab = "Forecast Year",
     main = paste("Forecast: Age-Specific Death Rate at age", x))

The hmd grabber functions, of which there are several in the R universe, typically do not check for errors returned from the server. Several HMD users report that hmd.mx() and friends broke after HMD did a server upgrade. In the server upgrade, the roles of www.mortality.org and mortality.org were switched. The server tried to tell the programs that the URL should be http://mortality.org instead of http://www.mortality.org.

Since this affected nearly all of the scripts for scraping R code, including one I wrote a long time ago, I changed the server back to the old way as recognizing www.mortality.org as the real server address, with mortality.org as a redirect. This may confuse a few browsers, but will save having code authors having to recode.

--Carl

I wish to use the code pop.sim without migration.
When trying this with the australian data provided by the addb package, I am getting the following error:

mort.sm <- smooth.demogdata(set.upperage(extract.years(australia,1950:2002),100))
fert.sm <- smooth.demogdata(extract.years(aus.fertility,1950:2002))
mort.fit <- coherentfdm(mort.sm)
fert.fit <- fdm(fert.sm)
mort.fcast <- forecast(mort.fit)
fert.fcast <- forecast(fert.fit)

aus.sim <- pop.sim(mort.fcast,fert.fcast,australia)
Error in pop.sim(mort.fcast, fert.fcast, australia) :
Inputs not fmforecast2 objects

I tried solving it by setting mig=NULL but then I am getting another error:
aus.sim <- pop.sim(mort.fcast,fert.fcast,mig=NULL,australia)
Error in pop.sim(mort.fcast, fert.fcast, mig = NULL, australia) :
Please ensure that mortality, migration and population data have the same age dimension

Is there any other way to solve this?
Thanks

Hello,

   I'm using the demography package to forecast Mauritian mortality rates. Mauritian data is not available in www.mortality.org. I have formatted my data in the same format as in Human mortality. But when I am trying to plot my data, it is giving a warning message : 

** total is my demogdata object

plot.demogdata(total)
Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

I'm having the same result for :

total.fit=fdm(total,order=4)
Error in if (max.age < max(ages)) data <- extract.ages(data, min(ages):max.age, :
missing value where TRUE/FALSE needed

If the data is 1x1 it works perfectly. However it does not work for data (5x1). It will be very grateful if you could help me. Thank you.

Hello,
With the new mortality.org web, can i still download data from HMD using this package?
Because, through my username and password using hmd.mx i get the following error Connection error at www.mortality.org. Please check username, password and country label. I am sure that username and password is correct. Can you help me?

Thank you.

Hi Rob,

I was wondering if it would be possible to download life expectancy at birth data from all countries in all years available. I do not know how to implement theses features in the library, but they would certainly be very welcomed !

hmd.e0( country= "ALL"
               , years= "ALL"
               , username= "myusername"
               , password= "mypswd")

The function hmd.mx does not connect with the HMD anymore.

When forecasting population with pop.sim in demography package, I find the dimnames in output confusing.
Here is my example:

library("demography")

mort.fo <- read.demogdata(file = "data/mort_no.txt", popfile = "data/pop_fo.txt", type = "mortality", label = "Føroyar", skip = 0)
mort.fo <- extract.years(mort.fo, 1985:2007)

fert.fo <- read.demogdata(file = "data/fert_rate.txt", popfile = "data/pop_fert.txt", type = "fertility", label = "Føroyar", skip = 0, lambda = 0.4)
fert.fo <- extract.years(fert.fo, 1985:2007)

mort.sm <- smooth.demogdata(extract.years(mort.fo, 1985:2006))
fert.sm <- smooth.demogdata(extract.years(fert.fo, 1985:2006))
mig.fo <- netmigration(mort.fo, fert.fo)

mort.fit <- coherentfdm(mort.sm)
fert.fit <- fdm(fert.sm)
mig.fit <- coherentfdm(mig.fo)

mort.fcast <- forecast(mort.fit, h = 40, level = 95)
fert.fcast <- forecast(fert.fit, h = 40, level = 95)
mig.fcast <- forecast(mig.fit, h = 40, level = 95)

set.seed(27052018) # dagfesting fyri fyrstu simulering
sim10000 <- pop.sim(mort = mort.fcast, fert = fert.fcast, mig = mig.fcast, firstyearpop = mort.fo, N = 10000, mfratio = 1.05, bootstrap = T)

Output

In this case the firstyearpop is medio 2007, and the first forecasted year have the dimname = 2007, it does not make sense to me, it should be 2008. Am I right or wrong?
I made a forecast 40 year a head. While comparing forecasted population with dimname = 2007, with actual population 2007, it does not fit at all, but fit perfectly with actual population 2008. Dimname 2016 fit also very well with actual population 2017.

What about changing pop.sim function to something like this:

at the very end of the pop.sim function, line number 153 and 154 it says:

dimnames(pop.m)[[2]] <- dimnames(pop.m)[[2]] <- firstyr +
(1:h) - 1
return(list(male = pop.m, female = pop.f))
}

Changed to:

dimnames(pop.m)[[2]] <- dimnames(pop.f)[[2]] <- firstyr +
(1:h)
return(list(male = pop.m, female = pop.f))
}

What do you think?

Thanks, cool work ! Package vignette would be extremely useful. The package seems to contain very interesting functionality - an online vignette would quickly give an overview what it really has to offer. At the moment the potential user has to browse the detailed function documentation which is missing visual examples.

when I try to load the package i get:

Load package: forecast
Error: package or namespace load failed for ‘demography’ in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]):
their is no package named ‘mclust’

I use R 3.5.3 on MacOs

I have no idea what to do.

Hi there

I'm using the demography package and I cannot run the lca function (it used to work in my old laptop, but since I have upgraded to a new one, it is not working anymore). Even if I try to run the first example in the help page
france.LC1 <- lca(fr.mort, adjust="e0")
I get the following error
Error in lca(fr.mort, adjust = "e0") : unused argument (adjust = "e0")

Any idea of what might be wrong? I'm using Windows 10 and the latest version of R and any contributed package

Best regards and thank you very much in advance

Javier

While using the command 'read.demogdata' and not setting any value for 'max.mx', I have noticed that the default value assigned to max.mx is 10 (This is not mentioned in the pdf file 'Package Demography V1.15'). Can someone confirm this for me?

the following code gives this error

hmd.mx(country = "JPN", username = "*********",password = "*******", label = "Japan") 

#####Error in hmd.mx(country = "JPN", username = "*********", password = "*******!",  : 
  Connection error at www.mortality.org. Please check username, password and country label.#####

Since Hmd have changed their interface website, I had to check-in with a new password but I still can't import data from there using hmd.mx

Reported at http://stats.stackexchange.com/q/54727/159

library(demography)
library(addb)
ausfert <- extract.years(aus.fert, years = 1921:2000)
ausf.sm <- smooth.demogdata(ausfert)
ausf.fit <- fdm(ausf.sm,order = 3, method = "classical")
plot(residuals(ausf.fit), type = "filled.contour")
plot(ausf.sm,transform=TRUE); lines(fitted(ausf.fit),col='black')

Hi,
I am having trouble compiling a pdf with Rmd while using the demography package for several weeks now.
The pdf compiles well but as soon as I add the hmd.mx() function it does not compile anymore. The R code with hmd.mx() works fine alone. I also tried creating a very simple Rmd with only that particular code and I get the same error. I also tried on a different computer. Here is the error message:

I was wondering if you have encountered a similar issue before and/or if you had some ideas for how to resolve it. I was also wondering if it could be a similar issue to the one in this post: colindouglas/retrosheet#1.

Thanks in advance for any help.

Rob,
The HMD website is now wrapped entirely in SSL so that passwords are no longer sent in the clear.

As a result, code using RCurl with a URL of http://www.mortality.org will get a 301 redirect returned, saying to use https://www.mortality.org

In my experiment with local edits of your hmd.mx() function, just changing http:// to https:// seemed to return the expected results.

Regards,
Carl Boe, HMD

running the code:

uk1<-hmd.mx("GBR_NP", login, password, "UK")

gives the following error message:

Error in read.table(con, header = TRUE, skip = 2, na.strings = ".", as.is = TRUE) :
more columns than column names

I have run hmd.mx function several times without encountering the issue before.

Thank you

A couple of refinements needed for e0 calculations - to do with minimum and maximum ages for life table calculations.

Usage

life.expectancy(data, series = names(data$rate)[1], years = data$year,
type = c("period", "cohort"), age = min(data$age),
max.age = min(100, max(data$age)))
flife.expectancy(data, series=NULL, years = data$year, type = c("period", "cohort"),
age = min(data$age), max.age = NULL, PI = FALSE, nsim = 500, ...)
e0(data, series = NULL, years = data$year, type = c("period", "cohort"),
max.age = NULL, PI = FALSE, nsim = 500, ...)

lifetable(data, series = names(data$rate)[1], years = data$year,
ages = data$age, max.age = min(100, max(data$age)),
type = c("period", "cohort"))

MINIMUM AGE is picked up from data in life.expectancy() and flife.expectancy() parameters but this is not used in e0()

e0
function (data, series = NULL, years = data$year, type = c("period",
"cohort"), max.age = NULL, PI = FALSE, nsim = 500, ...)
{
flife.expectancy(data, series = series, years = years, age = 0,
type = type, max.age = max.age, PI = PI, nsim = nsim,
...)
}

MAXIMUM AGE is not allowed to be greater than 100 by virtue of
max.age = min(100, max(data$age) in life.expectancy() and lifetable

If the user is happy with the smoothed rates at 110+, there is no reason not to allow the life table and life expectancy to use data to 110+. So specifying a max.age of 110 does not result in a forecast at ages 0-110+.

Grateful if specified minimum and maximum ages are carried through. Thanks

Prepare for release:

• git pull
• Check current CRAN check results
• Check if any deprecation processes should be advanced, as described in Gradual deprecation
• Polish NEWS
• devtools::build_readme()
• urlchecker::url_check()
• devtools::check(remote = TRUE, manual = TRUE)
• devtools::check_win_devel()
• rhub::check_for_cran()
• revdepcheck::revdep_check(num_workers = 4)
• Update cran-comments.md
• git push
• Draft blog post

Submit to CRAN:

• usethis::use_version('major')
• devtools::submit_cran()
• Approve email

Wait for CRAN...

• Accepted 🎉
• git push
• usethis::use_github_release()
• usethis::use_dev_version()
• git push
• Finish blog post
• Tweet
• Add link to blog post in pkgdown news menu

The code base was changed (commit 624f5f2) regarding dependency on library 'ftsa'. To DESCRIPTION was added 'Imports: ftsa (>= 4.8)' . I think that the intention was to remove 'import(ftsa)' from the NAMESPACE file, but that did not change.

library ftsa currently stands at version 3.8 on github, so the criterion in the DESCRIPTION file cannot be met and the demography package on github cannot be installed.

I am currently researching on the Lee-Carter(1992) model and hoping to make contribution to knowledge in the area of mortality forecasting. I went through the LCA function in R and I was able to produce a basic R summary codes for estimating the components of the Lee-Carter model (ax, bx and kt) .
However , I am wondering whether there is anything I am missing. My results for
"ax" seem to match with values from LCA function, however the values
for ,"bx" and "kt" have significant differences .
Could you kindly lend me some time to run the attached R codes. Is there anything I am missing? .
I would appreciate very well your input and suggestions.

REPLICATION_LCA FUNC.docx

When smoothing and passing ages beyond the range of the existing data, both the mortality rate and the population appear to be estimated with monotonic increasing splines.

I don't think having population should be monotonically increasing with age when extrapolating for older ages.

fr.2006 <- extract.years(fr.mort,2006)
fr.sm.2006 <- smooth.demogdata(fr.2006, age.grid = 0:115)
tail(fr.sm.2006$pop$total,10)
         2006
106 153.37000
107  71.66000
108  32.99000
109  13.71000
110   7.52000
111  14.29589
112  34.16178
113  67.11767
114 113.16356
115 172.29944

The country codes in the Human Mortality Database were changed. As a result, the country codes enumerated in the documentation of the function hmd.mx, as at 31 July 2012, are not up to date anymore.

When extracting data using hmd.mx I get the following error:

Error in function (type, msg, asError = TRUE) :
Unknown SSL protocol error in connection to www.mortality.org:443

I've used this hundreds of times, but probably not in last few months

The following causes problems.

fit.m <- fdm(fr.mort,'male')
fcast.m <- forecast(fit.m, h=100)
flife.expectancy(fcast.m, PI=TRUE,nsim=10, type="cohort", age=46, year=2007)

Leaving out year argument works.

Hi, I'm facing an issue when forecasting mortality data using Lee-Carter. This is the R code I used:

This is the re-adjusted kt output:

Thus, I would expect the forecast kt at year 2001 to be around -15. However, the generated value is -0.55:

After doing some digging, I found that the forecast function changed my re-adjusted kt value:

I'm not sure what went wrong since the ax and bx values are unaffected. Would appreciate on any advice to get the expected forecast kt value. Thanks!

Hello,
I'm trying to use the lca function for an analysis with information for the case of Switzerland. I start by downloading the data with the command
Mortality_info<-hmd.mx("CHE","username","password","CHE") (at this point I get a message saying something about NAs introduced for some convertion, but I get the message for other countries as well and I have no further issue.)

Later I just leave out older ages to avoid issues with nas/mortality errors through the command
Info_n<-extract.ages(Mortality_info,ages=0:90)

Then I try to run the function by using
res<-lca(Info_n,series="male",adjust="e0",years=1950:max(Info_n$year))
but I get the message:
Error in svd(clogrates) : infinite or missing values in 'x'

No matter if I do it for males or females or even if I change other parameters like years or ages, the function returns the same error though I see no missing values. I have tried to solve the problem unsuccessfully, and I do not really know what to do to be able to run the function properly. Do you have any ideas of why this is happening and how to solve it? It's really important.
Thanks in advance.

The following code snippet:

library(demography)

data <- hmd.mx("FRACNP",
               "[email protected]",
               "xxx",
               "France Civilian Population")

gives the following error message, although the username and password (in each case anonymised) and country label are correct:

Error in hmd.mx("FRACNP", "[email protected]", "xxx", "France Civilian Population") : 
  Connection error at www.mortality.org. Please check username, password and country label.

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] demography_1.20 forecast_8.4 

When there are lots of zeros in deaths, births and population numbers, the package can fail.

Need to: make existing functions more robust to zeros. (Test on Faroe Islands data).

Could add: interpolation using a bivariate smoother on age-time mortality surface.

Can I apply lca to abridged mortality rate? Age: 0, 1, 5, ...100 Years: 5 year steps. I tried it with 10 period data and following is what I get when I run the following codes

data <- read.demogdata(file="rate.txt", popfile="popfile.txt", type="mortality", label="NepalMale"
, max.mx = 10, skip = 0, popskip = 0, lambda=0, scale=1000)

nepalmale <- lca(data=data, series=names(data$rate)[1], years=data$year,
ages=data$age, max.age=100,
adjust = c("dt", "dxt", "e0", "none")[1], chooseperiod=FALSE,
minperiod=9, breakmethod=c("bai","bms"), scale = FALSE, restype = c("logrates", "rates", "deaths")[2],
interpolate = FALSE)

Error in matrix(mx, nrow = n, ncol = m) :
invalid 'ncol' value (too large or NA)

When loading the demography package, the ftsa package is required but it appears to require an installation of JAGS 3.00 > on the user machines. Is there any way to remove such dependency (maybe not attaching the package on load)

It would be nice to have a function (or a wrapper) that, given a fitted mortality projection model and well as its forecast, returns a vector of actual (and projected) mortality rates in order to obtain mortality tables to perform demographic and actuarial calculations.

Good day, Professor, I find it impossible to create demogdata object (as in example on your blog) if data is just population and not from hmd site (but has same structure).

Hello,
After the switch to the new site hmd.mx isn't working anymore, The error I'm getting is:

Error in function (type, msg, asError = TRUE) :
Unknown SSL protocol error in connection to www.mortality.org:443

HMDHFDplus is working fine so it's not a connection problem on my end.

Hello, through my username and password i was able to download mortality data from HMD using hmd.mx there is a week now that i get the following error Connection error at www.mortality.org. Please check username, password and country label. I am sure that username and password is correct. Can you help me?

Thank you.