smoothCluster: Cluster-level space-time smoothing models for mortality rates

Description

The function smoothCluster replace the previous function name fitINLA2 (before version 1.0.0).

Usage

smoothCluster(
  data,
  X = NULL,
  family = c("betabinomial", "binomial")[1],
  age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),
  age.n = c(1, 11, 12, 12, 12, 12),
  age.time.group = c(1, 2, 3, 3, 3, 3),
  age.strata.fixed.group = c(1, 2, 3, 4, 5, 6),
  time.model = c("rw1", "rw2", "ar1")[2],
  st.time.model = NULL,
  Amat,
  bias.adj = NULL,
  bias.adj.by = NULL,
  formula = NULL,
  year_label,
  type.st = 4,
  survey.effect = FALSE,
  linear.trend = TRUE,
  common.trend = FALSE,
  strata.time.effect = FALSE,
  hyper = "pc",
  pc.u = 1,
  pc.alpha = 0.01,
  pc.u.phi = 0.5,
  pc.alpha.phi = 2/3,
  pc.u.cor = 0.7,
  pc.alpha.cor = 0.9,
  pc.st.u = NA,
  pc.st.alpha = NA,
  pc.st.slope.u = NA,
  pc.st.slope.alpha = NA,
  overdisp.mean = 0,
  overdisp.prec = 0.4,
  options = list(config = TRUE),
  control.inla = list(strategy = "adaptive", int.strategy = "auto"),
  control.fixed = list(),
  verbose = FALSE,
  geo = NULL,
  rw = NULL,
  ar = NULL,
  st.rw = NULL,
  age.rw.group = NULL,
  ...
)
fitINLA2(
  data,
  X = NULL,
  family = c("betabinomial", "binomial")[1],
  age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),
  age.n = c(1, 11, 12, 12, 12, 12),
  age.time.group = c(1, 2, 3, 3, 3, 3),
  age.strata.fixed.group = c(1, 2, 3, 4, 5, 6),
  time.model = c("rw1", "rw2", "ar1")[2],
  st.time.model = NULL,
  Amat,
  bias.adj = NULL,
  bias.adj.by = NULL,
  formula = NULL,
  year_label,
  type.st = 4,
  survey.effect = FALSE,
  linear.trend = TRUE,
  common.trend = FALSE,
  strata.time.effect = FALSE,
  hyper = "pc",
  pc.u = 1,
  pc.alpha = 0.01,
  pc.u.phi = 0.5,
  pc.alpha.phi = 2/3,
  pc.u.cor = 0.7,
  pc.alpha.cor = 0.9,
  pc.st.u = NA,
  pc.st.alpha = NA,
  pc.st.slope.u = NA,
  pc.st.slope.alpha = NA,
  overdisp.mean = 0,
  overdisp.prec = 0.4,
  options = list(config = TRUE),
  control.inla = list(strategy = "adaptive", int.strategy = "auto"),
  control.fixed = list(),
  verbose = FALSE,
  geo = NULL,
  rw = NULL,
  ar = NULL,
  st.rw = NULL,
  age.rw.group = NULL,
  ...
)

Value

INLA model fit using the provided formula, country summary data, and geographic data

Arguments

data

count data of person-months with the following columns

cluster: cluster ID
years: time period
region: region of the cluster
strata: stratum of the cluster
age: age group corresponding to the row
total: total number of person-month in this age group, stratum, cluster, and period
Y: total number of deaths in this age group, stratum, cluster, and period

X

Covariate matrix. It must contain either a column with name "region", or a column with name "years", or both. The covariates must not have missing values for all regions (if varying in space) and all time periods (if varying in time). The rest of the columns are treated as covariates in the mean model.

family

family of the model. This can be either binomial (with logistic normal prior), betabiniomial.

age.groups

a character vector of age groups in increasing order.

age.n

number of months in each age groups in the same order.

age.time.group

vector indicating grouping of the ages groups in the temporal model. For example, if each age group is assigned a different temporal component, then set age.rw.group to c(1:length(age.groups)); if all age groups share the same random walk component, then set age.rw.group to a rep(1, length(age.groups)). The default for 6 age groups is c(1,2,3,3,3,3), which assigns a separate temporal trend to the first two groups and a common random walk for the rest of the age groups. The vector should contain values starting from 1. This argument replaces the previous age.rw.group argument.

age.strata.fixed.group

vector indicating grouping of the ages groups for different strata in the intercept. The default is c(1:length(age.groups)), which correspond to each age group within each stratum receives a separate intercept. If several age groups are specified to be the same value in this vector, the stratum specific deviation from the baseline is assumed to be the same for these age groups. For example, if age.strata.fixed.group = c(1, 2, 3, 3, 3, 3), then the intercept part of the linear predictor consists of 6 overall age-specific intercepts and 3 set of strata effects (where a base stratum is chosen internally), for age groups 1, 2, and the rest respectively. Notice that this argument does not control the linear temporal trends (which is also parameterized as fixed effect, but determined by the age.rw.group argument). The vector should contain values starting from 1.

More specific examples: (1) if each age group is assigned a different intercept, then set age.strata.fixed.group to c(1:length(age.groups)) (2) if all age groups share the same intercept, then set age.strata.fixed.group to a rep(1, length(age.groups)). The default for 6 age groups is the former. (3) If each temporal trend is associated with its own intercept, set it to be the same as age.rw.group.

time.model

Model for the main temporal trend, can be rw1, rw2, ar1, or NULL (for spatial-only smoothing). Default to be rw2. For ar1 main effect, a linear slope is also added with time scaled to be between -0.5 to 0.5, i.e., the slope coefficient represents the total change between the first year and the last year in the projection period on the logit scale.

st.time.model

Temporal component model for the interaction term, can be rw1, rw2, or ar1. Default to be the same as time.model unless specified otherwise. The default does not include region-specific random slopes. They can be added to the interaction term by specifying pc.st.slope.u and pc.st.slope.alpha.

Amat

Adjacency matrix for the regions

bias.adj

the ratio of unadjusted mortality rates or age-group-specific hazards to the true rates or hazards. It needs to be a data frame that can be merged to thee outcome, i.e., with the same column names for time periods (for national adjustment), or time periods and region (for subnational adjustment). The column specifying the adjustment ratio should be named "ratio".

bias.adj.by

vector of the column names specifying how to merge the bias adjustment to the count data. For example, if bias adjustment factor is provided in bias.adj for each region and time, then bias.adj.by should be `c("region", "time")`.

formula

INLA formula. See vignette for example of using customized formula.

year_label

string vector of year names

type.st

type for space-time interaction

survey.effect

logical indicator whether to include a survey fixed effect. If this is set to TRUE, there needs to be a column named 'survey' in the input data frame. In prediction, this effect term will be set to 0.

linear.trend

logical indicator whether a linear trend is added to the temporal main effect. If the temporal main effect is RW2, then it will be forced to FALSE. Default is TRUE.

common.trend

logical indicator whether all age groups and strata share the same linear trend in the temporal main effect.

strata.time.effect

logical indicator whether to include strata specific temporal trends.

hyper

Deprecated. which hyperpriors to use. Only supports PC prior ("pc").

pc.u

hyperparameter U for the PC prior on precisions.

pc.alpha

hyperparameter alpha for the PC prior on precisions.

pc.u.phi

hyperparameter U for the PC prior on the mixture probability phi in BYM2 model.

pc.alpha.phi

hyperparameter alpha for the PC prior on the mixture probability phi in BYM2 model.

pc.u.cor

hyperparameter U for the PC prior on the autocorrelation parameter in the AR prior, i.e. Prob(cor > pc.u.cor) = pc.alpha.cor.

pc.alpha.cor

hyperparameter alpha for the PC prior on the autocorrelation parameter in the AR prior.

pc.st.u

hyperparameter U for the PC prior on precisions for the interaction term.

pc.st.alpha

hyperparameter alpha for the PC prior on precisions for the interaction term.

pc.st.slope.u

hyperparameter U for the PC prior on precisions for the area-level random slope. If both pc.st.slope.u and pc.st.slope.alpha are not NA, an area-level random slope with iid prior will be added to the model. The parameterization of the random slope is so that Prob(|beta| > pc.st.slope.u) = pc.st.slope.alpha, where time covariate is rescaled to be -0.5 to 0.5, so that the random slope can be interpreted as the total deviation from the main trend from the first year to the last year to be projected, on the logit scale.

pc.st.slope.alpha

hyperparameter alpha for the PC prior on precisions for the area-level random slope. See above for the parameterization.

overdisp.mean

hyperparameter for the betabinomial likelihood. Mean of the over-dispersion parameter on the logit scale.

overdisp.prec

hyperparameter for the betabinomial likelihood. Precision of the over-dispersion parameter on the logit scale.

options

list of options to be passed to control.compute() in the inla() function.

control.inla

list of options to be passed to control.inla() in the inla() function. Default to the "adaptive" integration strategy.

control.fixed

list of options to be passed to control.fixed() in the inla() function.

verbose

logical indicator to print out detailed inla() intermediate steps.

geo

Deprecated. Spatial polygon file, legacy parameter from previous versions of the package.

rw

Deprecated. Take values 0, 1 or 2, indicating the order of random walk. If rw = 0, the autoregressive process is used instead of the random walk in the main trend. See the description of the argument ar for details.

ar

Deprecated. Order of the autoregressive component. If ar is specified to be positive integer, the random walk components will be replaced by AR(p) terms in the interaction part. The main temporal trend remains to be random walk of order rw unless rw = 0.

st.rw

Deprecated. Take values 1 or 2, indicating the order of random walk for the interaction term. If not specified, it will take the same order as the argument rw in the main effect. Notice that this argument is only used if ar is set to 0.

age.rw.group

Deprecated. Legacy parameter replaced by age.time.group.

...

arguments to be passed to the inla() function call.

Author

Zehang Richard Li

Examples

Run this code

if (FALSE) {
library(dplyr)
data(DemoData)
# Create dataset of counts
counts.all <- NULL
for(i in 1:length(DemoData)){
  counts <- getCounts(DemoData[[i]][, c("clustid", "time", "age", "died",
                                        "region", "strata")],
            variables = 'died', by = c("age", "clustid", "region", 
                                         "time", "strata"))
  counts <- counts %>% mutate(cluster = clustid, years = time, Y=died)
  counts$strata <- gsub(".*\\.","",counts$strata)
  counts$survey <- names(DemoData)[i] 
  counts.all <- rbind(counts.all, counts)
}

# fit cluster-level model on the periods
periods <- levels(DemoData[[1]]$time)
fit <- smoothCluster(data = counts.all, 
     Amat = DemoMap$Amat, 
     time.model = "rw2", 
     st.time.model = "rw1",
     strata.time.effect =  TRUE, 
     survey.effect = TRUE,
     family = "betabinomial",
     year_label = c(periods, "15-19"))
summary(fit)
est <- getSmoothed(fit, nsim = 1000)
plot(est$stratified, plot.CI=TRUE) + ggplot2::facet_wrap(~strata) 

# fit cluster-level space-time model with covariate
# notice without projected covariates, we use periods up to 10-14 only
# construct a random covariate matrix for illustration
periods <- levels(DemoData[[1]]$time)
X <- expand.grid(years = periods, 
       region = unique(counts.all$region))
X$X1 <- rnorm(dim(X)[1])
X$X2 <- rnorm(dim(X)[1])
fit.covariate <- smoothCluster(data = counts.all, 
   X = X,
     Amat = DemoMap$Amat, 
     time.model = "rw2", 
     st.time.model = "rw1",
     strata.time.effect =  TRUE, 
     survey.effect = TRUE,
     family = "betabinomial",
     year_label = c(periods))
est <- getSmoothed(fit.covariate, nsim = 1000)

# fit cluster-level model for one time point only
# i.e., space-only model
fit.sp <- smoothCluster(data = subset(counts.all, time == "10-14"), 
     Amat = DemoMap$Amat, 
     time.model = NULL, 
     survey.effect = TRUE,
     family = "betabinomial")
summary(fit.sp)
est <- getSmoothed(fit.sp, nsim = 1000)
plot(est$stratified, plot.CI = TRUE) + ggplot2::facet_wrap(~strata) 

# fit cluster-level model for one time point and covariate
# construct a random covariate matrix for illustration
X <- data.frame(region = unique(counts.all$region),
      X1 = c(1, 2, 2, 1), 
      X2 = c(1, 1, 1, 2))
fit.sp.covariate <- smoothCluster(data = subset(counts.all, time == "10-14"), 
     X = X, 
     Amat = DemoMap$Amat, 
     time.model = NULL, 
     survey.effect = TRUE,
     family = "betabinomial")
summary(fit.sp.covariate)
est <- getSmoothed(fit.sp.covariate, nsim = 1000)
}