Skip to content

Commit

Permalink
Ran styler across package
Browse files Browse the repository at this point in the history
  • Loading branch information
robjhyndman committed Dec 6, 2024
1 parent 8200edd commit 80e4d36
Show file tree
Hide file tree
Showing 16 changed files with 2,755 additions and 2,709 deletions.
18 changes: 9 additions & 9 deletions R/as.data.frame.demogdata.R
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,16 @@ as.data.frame.demogdata <- function(x, ...) {
# Size of matrices
nyears <- length(x$year)
nages <- length(x$age)
if(rates_included)
if (rates_included) {
groups <- names(x$rate)
else if(pop_included)
} else if (pop_included) {
groups <- names(x$pop)
else
} else {
groups <- NULL
outlist <- vector(length=nyears, mode="list")
}
outlist <- vector(length = nyears, mode = "list")
# Create data frame for rates
for(i in seq_along(groups)) {
for (i in seq_along(groups)) {
outlist[[i]] <- data.frame(
Year = rep(x$year, rep(nages, nyears)),
Age = rep(x$age, nyears),
Expand All @@ -40,7 +41,7 @@ as.data.frame.demogdata <- function(x, ...) {
out$Age <- as.integer(out$Age)
out$Year <- as.integer(out$Year)
# Assume Inf rates are due to 0/0
out$Rates[out$Rates==Inf] <- NA_real_
out$Rates[out$Rates == Inf] <- NA_real_
# Rename rates column
if (x$type == "mortality") {
colnames(out)[4] <- "Mortality"
Expand All @@ -62,10 +63,9 @@ as.data.frame.demogdata <- function(x, ...) {
}
}
# Reorganize
out <- out[order(out$Group, out$Year, out$Age),]
out <- out[order(out$Group, out$Year, out$Age), ]
rownames(out) <- NULL
return(out)
}

utils::globalVariables(c("Deaths","Births","Year","Age","Exposure","Group","Mortality","Fertility"))

utils::globalVariables(c("Deaths", "Births", "Year", "Age", "Exposure", "Group", "Mortality", "Fertility"))
5 changes: 2 additions & 3 deletions R/attach.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
.onAttach <- function(...)
{
msg <- paste("This is demography", utils::packageVersion("demography"),"\n")
.onAttach <- function(...) {
msg <- paste("This is demography", utils::packageVersion("demography"), "\n")
packageStartupMessage(msg)
}
187 changes: 94 additions & 93 deletions R/coherent.R
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

# COHERENT FORECASTING FOR MULTIPLE GROUPS

#' Coherent functional demographic model for grouped data
Expand Down Expand Up @@ -30,70 +29,74 @@
#' @seealso \code{\link{fdm}}, \code{\link{forecast.fdmpr}}
#'
#' @examples
#' fr.short <- extract.years(fr.sm,1950:2006)
#' fr.short <- extract.years(fr.sm, 1950:2006)
#' fr.fit <- coherentfdm(fr.short)
#' summary(fr.fit)
#' plot(fr.fit$product, components=3)
#' plot(fr.fit$product, components = 3)
#' @keywords models
#'
#'
#' @export
coherentfdm <- function(data, order1=6, order2=6, ...)
{
# Check if missing data

coherentfdm <- function(data, order1 = 6, order2 = 6, ...) {
# Check if missing data

# Use male and female if available
gps <- names(data$rate)
if(is.element("male",gps) & is.element("female",gps))
{
notneeded <- (1:length(gps))[-match(c("male","female"),gps)]
for(i in notneeded)
if (is.element("male", gps) & is.element("female", gps)) {
notneeded <- (1:length(gps))[-match(c("male", "female"), gps)]
for (i in notneeded) {
data$rate[[i]] <- data$pop[[i]] <- NULL
}
}

J <- length(data$rate)
rate.ratio <- fdm.ratio <- list()

# Construct ratio and product objects
is.mortality <- (data$type=="mortality")
is.mortality <- (data$type == "mortality")
y <- as.numeric(is.mortality)
pop.total <- 0
for (j in 1:J)
for (j in 1:J)
{
if(is.mortality)
y <- y * data$rate[[j]]^(1/J)
else
y <- y + data$rate[[j]]/J
if (is.mortality) {
y <- y * data$rate[[j]]^(1 / J)
} else {
y <- y + data$rate[[j]] / J
}
pop.total <- pop.total + data$pop[[j]]
}
rate.product <- demogdata(y, pop=pop.total, data$age, data$year, type=data$type, lambda=data$lambda,
label=data$label, name="product")
rate.product <- demogdata(y,
pop = pop.total, data$age, data$year, type = data$type, lambda = data$lambda,
label = data$label, name = "product"
)
for (j in 1:J)
{
if(is.mortality)
{
rate.ratio[[j]] <- demogdata(data$rate[[j]]/rate.product$rate[[1]], pop=pop.total, data$age,
data$year, type=data$type, label=data$label, name=names(data$rate)[[j]],lambda=data$lambda)
}
else
{
rate.ratio[[j]] <- demogdata(data$rate[[j]]-rate.product$rate[[1]], pop=pop.total, data$age,
data$year, type=data$type, label=data$label, name=names(data$rate)[[j]],lambda=data$lambda)
}
# Set infinite rates to missing
infrates <- rate.ratio[[j]]$rate[[1]] > 1e9
rate.ratio[[j]]$rate[[1]][infrates] <- NA
}

{
if (is.mortality) {
rate.ratio[[j]] <- demogdata(data$rate[[j]] / rate.product$rate[[1]],
pop = pop.total, data$age,
data$year, type = data$type, label = data$label, name = names(data$rate)[[j]], lambda = data$lambda
)
} else {
rate.ratio[[j]] <- demogdata(data$rate[[j]] - rate.product$rate[[1]],
pop = pop.total, data$age,
data$year, type = data$type, label = data$label, name = names(data$rate)[[j]], lambda = data$lambda
)
}
# Set infinite rates to missing
infrates <- rate.ratio[[j]]$rate[[1]] > 1e9
rate.ratio[[j]]$rate[[1]][infrates] <- NA
}

# GM model
fdm.mean <- fdm(rate.product, series="product", order=order1, ...)
fdm.mean <- fdm(rate.product, series = "product", order = order1, ...)

# Ratio model
for (j in 1:J)
fdm.ratio[[j]] <- fdm(rate.ratio[[j]], series=names(data$rate)[j], order=order2, ...)
for (j in 1:J) {
fdm.ratio[[j]] <- fdm(rate.ratio[[j]], series = names(data$rate)[j], order = order2, ...)
}
names(fdm.ratio) <- names(data$rate)
return(structure(list(product=fdm.mean, ratio=fdm.ratio), class="fdmpr"))

return(structure(list(product = fdm.mean, ratio = fdm.ratio), class = "fdmpr"))
}


Expand All @@ -118,91 +121,89 @@ coherentfdm <- function(data, order1=6, order2=6, ...)
#' finally a list of forecasts from each of the ratio models.
#' @author Rob J Hyndman
#' @seealso \code{\link{coherentfdm}}, \code{\link{forecast.fdm}}.
#' @examples fr.short <- extract.years(fr.sm,1950:2006)
#' @examples fr.short <- extract.years(fr.sm, 1950:2006)
#' fr.fit <- coherentfdm(fr.short)
#' fr.fcast <- forecast(fr.fit)
#' plot(fr.fcast$male)
#' plot(fr.fcast$ratio$male, plot.type='component', components=3)
#' plot(fr.fcast$ratio$male, plot.type = "component", components = 3)
#' models(fr.fcast)
#'
#'
#' @keywords models
#' @export
forecast.fdmpr <- function(object, h=50, level=80, K=100, drange=c(0.0,0.5), ...)
{
forecast.fdmpr <- function(object, h = 50, level = 80, K = 100, drange = c(0.0, 0.5), ...) {
fcast.ratio <- fc <- totalvar.r <- list()
J <- length(object$ratio)
ny <- length(object$ratio[[1]]$year)
K <- min(K,ny)
K <- min(K, ny)

# GM model
fcast.mean <- forecast(object$product, method="arima", h=h, level=level, ...)
fcast.mean <- forecast(object$product, method = "arima", h = h, level = level, ...)
# Make sure first coefficient is not I(1) with drift.
#mod <- auto.arima(object$product$coeff[,2],d=2)
#fcast.mean$coeff[[2]] <- forecast(mod, h=h, level=level, ...)
#fcast.mean <- update(fcast.mean)
# mod <- auto.arima(object$product$coeff[,2],d=2)
# fcast.mean$coeff[[2]] <- forecast(mod, h=h, level=level, ...)
# fcast.mean <- update(fcast.mean)

# Obtain forecasts for each group
is.mortality <- (object$product$type=="mortality")
y <- as.numeric(is.mortality) #=1 for mortality and 0 for migration
for (j in 1:J)
is.mortality <- (object$product$type == "mortality")
y <- as.numeric(is.mortality) # =1 for mortality and 0 for migration
for (j in 1:J)
{
# Use all available data other than last K years
# As ARFIMA can't handle missing values
object$ratio[[j]]$weights <- 0*object$ratio[[j]]$weights +1
if(K < ny)
object$ratio[[j]]$weights[1:(ny-K)] <- 0
fcast.ratio[[j]] <- forecast(object$ratio[[j]], h=h, level=level, method="arfima", estim="mle", drange=drange,...)
object$ratio[[j]]$weights <- 0 * object$ratio[[j]]$weights + 1
if (K < ny) {
object$ratio[[j]]$weights[1:(ny - K)] <- 0
}
fcast.ratio[[j]] <- forecast(object$ratio[[j]], h = h, level = level, method = "arfima", estim = "mle", drange = drange, ...)
fc[[j]] <- fcast.mean
if(is.mortality)
if (is.mortality) {
fc[[j]]$rate[[1]] <- fcast.mean$rate$product * fcast.ratio[[j]]$rate[[1]]
else
} else {
fc[[j]]$rate[[1]] <- fcast.mean$rate$product + fcast.ratio[[j]]$rate[[1]]
}
names(fc[[j]]$rate)[1] <- names(object$ratio)[j]
fc[[j]]$coeff <- fc[[j]]$coeff.error <- fc[[j]]$call <- fc[[j]]$var <- NULL
if(is.mortality)
if (is.mortality) {
y <- y * fc[[j]]$rate[[1]]
else
} else {
y <- y + fc[[j]]$rate[[1]]
}
}

# Adjust forecasts so they multiply appropriately.
if(is.mortality)
{
y <- y^(1/J)/fcast.mean$rate$product
for(j in 1:J)
fc[[j]]$rate[[1]] <- fc[[j]]$rate[[1]]/y
}
else
{
y <- y/J - fcast.mean$rate$product
for(j in 1:J)
fc[[j]]$rate[[1]] <- fc[[j]]$rate[[1]]-y
}
if (is.mortality) {
y <- y^(1 / J) / fcast.mean$rate$product
for (j in 1:J) {
fc[[j]]$rate[[1]] <- fc[[j]]$rate[[1]] / y
}
} else {
y <- y / J - fcast.mean$rate$product
for (j in 1:J) {
fc[[j]]$rate[[1]] <- fc[[j]]$rate[[1]] - y
}
}
# Variance of forecasts
qconf <- 2 * stats::qnorm(0.5 + fcast.mean$coeff[[1]]$level/200)
for (j in 1:J)
qconf <- 2 * stats::qnorm(0.5 + fcast.mean$coeff[[1]]$level / 200)
for (j in 1:J)
{
vartotal <- fcast.mean$var$total + fcast.ratio[[j]]$var$total
tmp <- qconf * sqrt(vartotal)
fc[[j]]$rate$lower <- InvBoxCox(BoxCox(fc[[j]]$rate[[1]],object$product$lambda) - tmp, object$product$lambda)
fc[[j]]$rate$upper <- InvBoxCox(BoxCox(fc[[j]]$rate[[1]],object$product$lambda) + tmp, object$product$lambda)
if(is.mortality)
{
fc[[j]]$model[[4]] <- BoxCox(InvBoxCox(object$product[[4]],object$product$lambda) *
InvBoxCox(object$ratio[[j]][[4]], object$product$lambda),object$product$lambda)
}
else
{
fc[[j]]$model[[4]] <- BoxCox(InvBoxCox(object$product[[4]],object$product$lambda) +
InvBoxCox(object$ratio[[j]][[4]], object$product$lambda),object$product$lambda)
fc[[j]]$rate$lower <- InvBoxCox(BoxCox(fc[[j]]$rate[[1]], object$product$lambda) - tmp, object$product$lambda)
fc[[j]]$rate$upper <- InvBoxCox(BoxCox(fc[[j]]$rate[[1]], object$product$lambda) + tmp, object$product$lambda)
if (is.mortality) {
fc[[j]]$model[[4]] <- BoxCox(InvBoxCox(object$product[[4]], object$product$lambda) *
InvBoxCox(object$ratio[[j]][[4]], object$product$lambda), object$product$lambda)
} else {
fc[[j]]$model[[4]] <- BoxCox(InvBoxCox(object$product[[4]], object$product$lambda) +
InvBoxCox(object$ratio[[j]][[4]], object$product$lambda), object$product$lambda)
}
names(fc[[j]]$model)[4] <- names(object$ratio)[j]
fc[[j]]$coeff <- list(list(level= fcast.mean$coeff[[1]]$level))
fc[[j]]$coeff <- list(list(level = fcast.mean$coeff[[1]]$level))
}

names(fc) <- names(fcast.ratio) <- names(object$ratio)
fc$product <- fcast.mean
fc$ratio <- fcast.ratio
return(structure(fc, class="fmforecast2"))

return(structure(fc, class = "fmforecast2"))
}
Loading

0 comments on commit 80e4d36

Please sign in to comment.