Package 'SUMMER'

Description

SUMMER provides methods for spatial and spatio-temporal smoothing of demographic and health indicators using survey data, with particular focus on estimating and projecting under-five mortality rates.

Details

For details on the model implemented in this package, see Mercer et al. (2015) <doi:10.1214/15-AOAS872> and Li et al. (2019) <doi:10.1371/journal.pone.0210645>.

The development version of the package will be maintained on https://github.com/richardli/SUMMER.

Helper function of pixelPopToArea

Description

Aggregates population from the pixel level to the level of the area of interest.

Usage

aggPixelPreds(
  Zg,
  Ng,
  areas,
  urban = targetPopMat$urban,
  targetPopMat = NULL,
  useDensity = FALSE,
  stratifyByUrban = TRUE,
  normalize = useDensity
)

Arguments

Description

Takes simulated populations and aggregates them to the specified areal level. Also calculates the aggregated risk and prevalence.

Usage

pixelPopToArea(
  pixelLevelPop,
  eaSamples,
  areas,
  stratifyByUrban = TRUE,
  targetPopMat = NULL,
  doFineScaleRisk = !is.null(pixelLevelPop$fineScaleRisk$p),
  doSmoothRisk = !is.null(pixelLevelPop$smoothRisk$p)
)

areaPopToArea(
  areaLevelPop,
  areasFrom,
  areasTo,
  stratifyByUrban = TRUE,
  doFineScaleRisk = !is.null(areaLevelPop$aggregationResults$pFineScaleRisk),
  doSmoothRisk = !is.null(areaLevelPop$aggregationResults$pSmoothRisk)
)

Arguments

Value

A list containing elements 'fineScalePrevalence' and 'fineScaleRisk'. Each of these are in turn lists with aggregated prevalence and risk for the area of interest, containg the following elements, were paranethesis indicate the elements for the fineScaleRisk model rather than fineScalePrevalence:

Functions

• pixelPopToArea(): Aggregate from pixel to areal level

• areaPopToArea(): Aggregate areal populations to another areal level

Author(s)

John Paige

References

In Preparation

See Also

areaPopToArea

Examples

## Not run: 
##### Now we make a model for the risk. We will use an SPDE model with these 
##### parameters for the linear predictor on the logist scale, which are chosen 
##### to be of practical interest:
beta0=-2.9 # intercept
gamma=-1 # urban effect
rho=(1/3)^2 # spatial variance
effRange = 400 # effective spatial range in km
sigmaEpsilon=sqrt(1/2.5) # cluster (nugget) effect standard deviation

# simulate the population! Note that this produces multiple dense 
# nEA x nsim and nIntegrationPoint x nsim matrices. In the future 
# sparse matrices will and chunk by chunk computations may be incorporated.
simPop = simPopSPDE(nsim=1, easpa=easpaKenyaNeonatal, 
                    popMat=popMatKenya, targetPopMat=popMatKenyaNeonatal, 
                    poppsub=poppsubKenya, spdeMesh=kenyaMesh, 
                    margVar=rho, sigmaEpsilonSq=sigmaEpsilon^2, 
                    gamma=gamma, effRange=effRange, beta0=beta0, 
                    seed=123, inla.seed=12, nHHSampled=25, 
                    stratifyByUrban=TRUE, subareaLevel=TRUE, 
                    doFineScaleRisk=TRUE, 
                    min1PerSubarea=TRUE)

pixelPop = simPop$pixelPop
subareaPop = pixelPopToArea(pixelLevelPop=pixelPop, eaSamples=pixelPop$eaSamples, 
  areas=popMatKenya$subarea, stratifyByUrban=TRUE, 
  targetPopMat=popMatKenyaNeonatal, doFineScaleRisk=TRUE)

# get areas associated with each subarea for aggregation
tempAreasFrom = popMatKenya$subarea
tempAreasTo = popMatKenya$area
areasFrom = sort(unique(tempAreasFrom))
areasToI = match(areasFrom, tempAreasFrom)
areasTo = tempAreasTo[areasToI]

# do the aggregation from subareas to areas
outAreaLevel = areaPopToArea(areaLevelPop=subareaPop, 
  areasFrom=areasFrom, areasTo=areasTo, 
  stratifyByUrban=TRUE, doFineScaleRisk=TRUE)

## End(Not run)

Description

Aggregate estimators from different surveys.

Usage

aggregateSurvey(data)

Arguments

Value

Estimators aggregated across surveys.

Author(s)

Zehang Richard Li

Examples

## Not run: 
data(DemoData)
data(DemoMap)
years <- levels(DemoData[[1]]$time)

# obtain direct estimates
data <- getDirectList(births = DemoData, 
years = years, 
regionVar = "region", timeVar = "time", 
clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", 
geo.recode = NULL)

# obtain maps
geo <- DemoMap$geo
mat <- DemoMap$Amat

# Simulate hyper priors
priors <- simhyper(R = 2, nsamp = 1e+05, nsamp.check = 5000, Amat = mat, only.iid = TRUE)

# combine data from multiple surveys
data <- aggregateSurvey(data)
utils::head(data)


## End(Not run)

Description

Benchmark posterior draws to national estimates

Usage

Benchmark(
  fitted,
  national,
  estVar,
  sdVar,
  timeVar = NULL,
  weight.region = NULL,
  method = c("MH", "Rejection")[2]
)

Arguments

Value

Benchmarked object in S3 class SUMMERproj or SUMMERprojlist in the same format as the input object fitted.

Author(s)

Taylor Okonek, Zehang Richard Li

Examples

## Not run: 
##  ------------------------------------------ ##
##     Benchmarking with smoothCluster output
##  ------------------------------------------ ##

data(DemoData)
# fit unstratified cluster-level model
counts.all <- NULL
for(i in 1:length(DemoData)){
vars <- c("clustid", "region", "time", "age")
counts <- getCounts(DemoData[[i]][, c(vars, "died")], 
						variables = 'died',
	 			    by = vars, drop=TRUE)
counts$cluster <- counts$clustid
counts$years <- counts$time
counts$Y <- counts$died
counts$survey <- names(DemoData)[i]	
counts.all <- rbind(counts.all, counts)
}
periods <- c("85-89", "90-94", "95-99", "00-04", "05-09", "10-14", "15-19")
fit.bb  <- smoothCluster(data = counts.all, Amat = DemoMap$Amat, 
				family = "betabinomial",
				year_label = periods, 
				survey.effect = TRUE)
est.bb <- getSmoothed(fit.bb, nsim = 1e4, CI = 0.95, save.draws=TRUE)

# construct a simple population weight data frame with equal weights
weight.region <- expand.grid(region = unique(counts.all$region), 
						 years = periods)
weight.region$proportion <- 0.25

# construct a simple national estimates
national <- data.frame(years = periods, 
				   est = seq(0.27, 0.1, length = 7), 
				   sd = runif(7, 0.01, 0.03))

 # benchmarking
est.bb.bench <- Benchmark(est.bb, national, weight.region = weight.region, 
						estVar = "est", sdVar = "sd", timeVar = "years")

# Sanity check: Benchmarking comparison
compare <- national
compare$before <- NA
compare$after <- NA
for(i in 1:dim(compare)[1]){
	weights <- subset(weight.region, years == national$years[i])
	sub <- subset(est.bb$overall, years == national$years[i])
	sub <- merge(sub, weights)
	sub.bench <- subset(est.bb.bench$overall, years == national$years[i])
	sub.bench <- merge(sub.bench, weights)
	compare$before[i] <- sum(sub$proportion * sub$median)
	compare$after[i] <- sum(sub.bench$proportion * sub.bench$median)
}
plot(compare$est, compare$after, col = 2, pch = 10,
		 xlim = range(c(compare$est, compare$before, compare$after)),
		 ylim = range(c(compare$est, compare$before, compare$after)),
		 xlab = "External national estimates", 
		 ylab = "Weighted posterior median after benchmarking",
    main = "Sanity check: weighted average of area medians")
points(compare$est, compare$before)
abline(c(0, 1))
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))

#  construct a simple national estimates
national <- data.frame(years = periods, 
					   est = seq(0.22, 0.1, length = 7), 
					   sd = runif(7, 0.01, 0.03))
# national does not need to have all years
national_sub <- national[1:3,]

# benchmarking
est.bb.bench <- Benchmark(est.bb, national_sub, 
						weight.region = weight.region, 
						estVar = "est", sdVar = "sd", timeVar = "years")

# Sanity check: only benchmarked for three periods
compare <- national
compare$before <- NA
compare$after <- NA
for(i in 1:dim(compare)[1]){
	weights <- subset(weight.region, years == national$years[i])
	sub <- subset(est.bb$overall, years == national$years[i])
	sub <- merge(sub, weights)
	sub.bench <- subset(est.bb.bench$overall, years == national$years[i])
	sub.bench <- merge(sub.bench, weights)
	compare$before[i] <- sum(sub$proportion * sub$median)
	compare$after[i] <- sum(sub.bench$proportion * sub.bench$median)
}
plot(compare$est, compare$after, col = 2, pch = 10,
		 xlim = range(c(compare$est, compare$before, compare$after)),
		 ylim = range(c(compare$est, compare$before, compare$after)),
		 xlab = "External national estimates", 
		 ylab = "Weighted posterior median after benchmarking",
    main = "Sanity check: weighted average of area medians")
points(compare$est, compare$before)
abline(c(0, 1))
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))

#  Another extreme benchmarking example, where almost all weights in central region
weight.region$proportion <- 0.01
weight.region$proportion[weight.region$region == "central"] <- 0.97
# benchmarking
est.bb.bench <- Benchmark(est.bb, national, weight.region = weight.region, 
					estVar = "est", sdVar = "sd", timeVar = "years")
# It can be seen the central region are pulled to the national benchmark
plot(national$est, 
	 subset(est.bb.bench$overall, region == "central")$mean,
	 col = 2, pch = 10, xlab = "External national estimates", 
	 ylab = "Central region estimates") 
points(national$est, 
	 subset(est.bb$overall, region == "central")$mean) 
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10),  col = c(1, 2))
abline(c(0, 1))

# Example with the MH method
# Benchmarking with MH should be applied when customized priors are 
#  specified for fixed effects when fitting the model
fit.bb.new  <- smoothCluster(data = counts.all, Amat = DemoMap$Amat, 
				family = "betabinomial",
				year_label = periods, 
				survey.effect = TRUE, 
				control.fixed = list(
					mean=list(`age.intercept0:1`=-4, 
						       `age.intercept1-11:1`=-5,
						       `age.intercept12-23:1`=-8,
						       `age.intercept24-35:1`=-9,
						       `age.intercept36-47:1`=-10,
						       `age.intercept48-59:1`=-11), 
					prec=list(`age.intercept0:1`=10, 
						       `age.intercept1-11:1`=10,
						       `age.intercept12-23:1`=10,
						       `age.intercept24-35:1`=10,
						       `age.intercept36-47:1`=10,
						       `age.intercept48-59:1`=10)))
est.bb.new <- getSmoothed(fit.bb.new, nsim = 10000, CI = 0.95, save.draws=TRUE)

#  construct a simple national estimates
national <- data.frame(years = periods, 
					   est = seq(0.22, 0.1, length = 7), 
					   sd = runif(7, 0.01, 0.03))
weight.region <- expand.grid(region = unique(counts.all$region), 
						 years = periods)
weight.region$proportion <- 0.25					   
est.bb.bench.MH <- Benchmark(est.bb.new, national, 
	weight.region = weight.region, 
	estVar = "est", sdVar = "sd", timeVar = "years",
	method = "MH")

compare <- national
compare$before <- NA
compare$after <- NA
for(i in 1:dim(compare)[1]){
	weights <- subset(weight.region, years == national$years[i])
	sub <- subset(est.bb.new$overall, years == national$years[i])
	sub <- merge(sub, weights)
	sub.bench <- subset(est.bb.bench.MH$overall, years == national$years[i])
	sub.bench <- merge(sub.bench, weights)
	compare$before[i] <- sum(sub$proportion * sub$median)
	compare$after[i] <- sum(sub.bench$proportion * sub.bench$median)
}
plot(compare$est, compare$after, col = 2, pch = 10,
		 xlim = range(c(compare$est, compare$before, compare$after)),
		 ylim = range(c(compare$est, compare$before, compare$after)),
		 xlab = "External national estimates", 
		 ylab = "Weighted posterior median after benchmarking",
    main = "Sanity check: weighted average of area medians")
points(compare$est, compare$before)
abline(c(0, 1))
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))

##  ------------------------------------------ ##
##     Benchmarking with smoothDirect output
##  ------------------------------------------ ##
years <- levels(DemoData[[1]]$time)
# obtain direct estimates
data_multi <- getDirectList(births = DemoData, years = years,
                        regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
                        ageVar = "age", weightsVar = "weights", geo.recode = NULL)
data <- aggregateSurvey(data_multi)
#  subnational model
years.all <- c(years, "15-19")
fit2 <- smoothDirect(data = data, Amat = DemoMap$Amat,
                 year_label = years.all, year_range = c(1985, 2019),
                 time.model = "rw2", m = 5, type.st = 4)
out2a <- getSmoothed(fit2, joint = TRUE, nsim = 1e5, save.draws = TRUE)

##
## Benchmarking for yearly estimates
##
weight.region <- expand.grid(region = unique(data$region[data$region != "All"]),
                             years = 1985:2019)
weight.region$proportion <- 0.25
# construct a simple national estimates
national <- data.frame(years = 1985:2019,
                       est = seq(0.25, 0.15, length = 35),
                       sd = runif(35, 0.03, 0.05))
# Benchmarking to national estimates on the yearly scale
out2b <- Benchmark(out2a, national, weight.region = weight.region,
                          estVar = "est", sdVar = "sd", timeVar = "years")
plot(out2a$overall)  
plot(out2b$overall) 

# combine the point estimate and compare with the benchmark values
national.est <- aggregate(mean ~ years, 
   data = out2a$overall[out2a$overall$is.yearly, ], FUN = mean)
national.est.bench <- aggregate(mean ~ years, 
   data = out2b$overall[out2b$overall$is.yearly, ], FUN = mean)

plot(national$est, national.est$mean,  
		 xlim = range(c(national$est, national.est$mean, national.est.bench$mean)),
		 ylim = range(c(national$est, national.est$mean, national.est.bench$mean)),
		 xlab = "External national estimates", 
		 ylab = "Weighted posterior median after benchmarking",
    main = "Sanity check: weighted average of area means")
points(national$est, national.est.bench$mean, col = 2, pch = 10)
abline(c(0, 1))
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))



##
## Benchmarking for period estimates
##
weight.region <- expand.grid(region = unique(data$region[data$region != "All"]),
                             years = years.all)
weight.region$proportion <- 0.25
# construct a simple national estimates
national <- data.frame(years = years.all,
                       est = seq(0.25, 0.15, len = 7),
                       sd = runif(7, 0.01, 0.03))
# Benchmarking to national estimates on the period scale
out2c <- Benchmark(out2a, national, weight.region = weight.region,
                          estVar = "est", sdVar = "sd", timeVar = "years")
plot(out2a$overall)
plot(out2c$overall)

# combine the point estimate and compare with the benchmark values
national.est <- aggregate(mean ~ years, 
			data = out2a$overall[!out2a$overall$is.yearly, ], FUN = mean)
national.est.bench <- aggregate(mean ~ years, 
			data = out2c$overall[!out2b$overall$is.yearly, ], FUN = mean)

plot(national$est, national.est$mean,  
		 xlim = range(c(national$est, national.est$mean, national.est.bench$mean)),
		 ylim = range(c(national$est, national.est$mean, national.est.bench$mean)),
		 xlab = "External national estimates", 
		 ylab = "Weighted posterior median after benchmarking",
    main = "Sanity check: weighted average of area means")
points(national$est, national.est.bench$mean, col = 2, pch = 10)
abline(c(0, 1))
legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))


 
## End(Not run)

Description

The Behavioral Risk Factor Surveillance System (BRFSS) is an annual telephone health survey conducted by the Centers for Disease Control and Prevention (CDC) that tracks health conditions and risk behaviors in the United States and its territories since 1984. This BRFSS dataset contains 16124 observations. The 'diab2' variable is the binary indicator of Type II diabetes, 'strata' is the strata indicator and 'rwt_llcp' is the final design weight. Records with missing HRA code or diabetes status are removed from this dataset. See https://www.cdc.gov/brfss/annual_data/2013/pdf/Weighting_Data.pdf for more details of the weighting procedure.

Usage

data(BRFSS)

Format

A data.frame of 26 variables.

References

Washington State Department of Health, Center for Health Statistics. Behavioral Risk Factor Surveillance System, supported in part by the Centers for Disease Control and Prevention. Corporative Agreement U58/DP006066-01 (2015).

Description

Calibrate/normalize the point level totals so their sum matches the regional totals. Technically, the totals can be at any level smaller than the region level specified.

Usage

calibrateByRegion(pointTotals, pointRegions, regions, regionTotals)

Arguments

Value

A vector of same length as pointTotals and pointRegions containing the calibrated/normalized point totals that sum to the correct regional totals

Vector of updated point level totals, calibrated to match region totals

Author(s)

John Paige

Examples

pointTotals = c(1, 1, 1, 2)
pointRegions = c("a", "a", "b", "b")
regionTotals = c(10, 20)
regions = c("a", "b")
calibrateByRegion(pointTotals, pointRegions, regions, regionTotals)

Description

Map region names to a common set.

Usage

ChangeRegion(data, Bmat, regionVar = "region")

Arguments

Value

Data after changing region names

Author(s)

Zehang Richard Li

Examples

# Construct a small test data
testdata <- data.frame(region = c("north", "south", "east",
 "south", "east"), index = c(1:5))

# Construct a changing rule: combining south and east
Bmat <- matrix(c(1, 0, 0, 0, 1, 1), 3, 2)
colnames(Bmat) <- c("north", "south and east")
rownames(Bmat) <- c("north", "south", "east")
print(Bmat)

# New data after transformation
test <- ChangeRegion(testdata, Bmat, "region")
print(test)

Description

Plot heatmap comparing pairwise posterior exceedence probabilities for svysae object

Usage

compareEstimates(x, posterior.sample = NULL, title = NULL, return.plot = FALSE)

Arguments

Value

ggplot containing heat map of pairwise comparisons

Examples

## Not run: 
data(DemoData2)
data(DemoMap2)
library(survey)
des0 <- svydesign(ids = ~clustid+id, strata = ~strata,
                  weights = ~weights, data = DemoData2, nest = TRUE)
Xmat <- aggregate(age~region, data = DemoData2, FUN = mean)

cts.res <- smoothArea(tobacco.use ~ 1,
                      domain = ~region,
                      design = des0,
                      adj.mat = DemoMap2$Amat, 
                      pc.u = 1,
                      pc.alpha = 0.01,
                      pc.u.phi = 0.5,
                      pc.alpha.phi = 2/3,
                      return.samples = TRUE)
compareEstimates(cts.res)

## End(Not run)

Description

A small simulated dataset with 4 regions and 5 survey years. This does not represent any real country's data and are based on a subset of the model dataset provided by DHS.

Usage

data(DemoData)

Format

A list of with five components, named by survey year.

Source

https://dhsprogram.com/data/model-datasets.cfm

Description

A small fake dataset with 8 regions and two response variables: age and tobacco.use. This does not represent any real country's data and are based on a subset of the model dataset provided by DHS.

Usage

data(DemoData2)

Format

A data.frame of 7 variables.

Source

https://dhsprogram.com/data/model-datasets.cfm

Description

Shapefiles are from 1995 Uganda Admin 1 regions provided by DHS, but the data do not represent real information about any country.

Usage

data(DemoMap)

Format

An object of class list of length 2.

Details

• geo. Geographic map files

• Amat. Adjacency matrix for regions

Source

https://spatialdata.dhsprogram.com/boundaries/#view=table&countryId=UG

Description

Shapefiles are from 2014 Kenya Admin 1 regions provided by DHS.

Usage

data(DemoMap2)

Format

An object of class list of length 2.

Details

• geo Geographic map files

• Amat Adjacency matrix for regions

Source

https://spatialdata.dhsprogram.com/boundaries/#view=table&countryId=KE

Description

Expit transformation

Usage

expit(x)

Arguments

Value

expit of x

Examples

x <- .5
expit(x)

Description

Adjust direct estimates and their associated variances

Usage

getAdjusted(
  data,
  ratio,
  time = "years",
  region = "region",
  est = "mean",
  logit = "logit.est",
  logit.var = "var.est",
  logit.prec = "logit.prec",
  logit.lower = "lower",
  logit.upper = "upper",
  prob.lower = NULL,
  prob.upper = NULL,
  adj = "ratio",
  verbose = FALSE,
  lower = NULL,
  upper = NULL
)

Arguments

Value

adjusted dataset of the same columns.

Author(s)

Zehang Richard Li

Examples

## Not run: 
years <- levels(DemoData[[1]]$time)

# obtain direct estimates
data <- getDirectList(births = DemoData, 
years = years,  
regionVar = "region", timeVar = "time", 
clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", 
geo.recode = NULL)
# obtain direct estimates
data_multi <- getDirectList(births = DemoData, years = years,
  regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
  ageVar = "age", weightsVar = "weights", geo.recode = NULL)
data <- aggregateSurvey(data_multi)

# randomly simulate adjustment factor
adj <- expand.grid(region = unique(data$region), years = years)
adj$ratio <- runif(dim(adj)[1], min = 0.5, max = 0.8)
data.adj <- getAdjusted(data = data, ratio = adj)
 
## End(Not run)

Description

Extract adjacency matrix from the map

Usage

getAmat(geo, names)

Arguments

Value

Spatial djacency matrix.

Author(s)

Zehang Richard Li

Examples

## Not run: 
data(DemoMap) 
mat <- getAmat(geo = DemoMap$geo, names = DemoMap$geo$REGNAME)
mat
DemoMap$Amat

## End(Not run)

Description

For any points not in an area, they are assigned the nearest area using fields::fields.rdist.near or fields::rdist depending on the number of points and the maximum memory in bytes with a warning.

Usage

getAreaName(
  pts,
  shapefile,
  areaNameVar = "NAME_1",
  delta = 0.05,
  mean.neighbor = 50,
  maxBytes = 3 * 2^30
)

Arguments

Details

delta and mean.neighbor arguments only used when some points are not in areas, perhaps due to inconsistencies in shapefiles.

Value

A list of area IDs, area names, whether or not points are in multiple areas, and whether or not points are in no areas and assigned to the nearest one.

Author(s)

John Paige

See Also

projKenya, fields.rdist.near

Examples

## Not run: 
# download Kenya GADM shapefiles from SUMMERdata github repository
githubURL <- "https://github.com/paigejo/SUMMERdata/blob/main/data/kenyaMaps.rda?raw=true"
download.file(githubURL,"kenyaMaps.rda")

# load it in
load("kenyaMaps.rda")

# use the shapefile data to see what Admin1 and 2 areas the 
# points (0, 37) and (0.5, 38) are in
# (these are longitude/latitude coordinates)
pts = cbind(c(37, 38), c(0, .5))
head(slot(adm1, "data"))
admin1Areas = getAreaName(pts, adm1, "NAME_1")
admin2Areas = getAreaName(pts, adm2, "NAME_2")

## End(Not run)

Description

Reformat full birth records into person-month format

Usage

getBirths(
  filepath = NULL,
  data = NULL,
  surveyyear = NA,
  variables = c("caseid", "v001", "v002", "v004", "v005", "v021", "v022", "v023", "v024",
    "v025", "v139", "bidx"),
  strata = c("v024", "v025"),
  dob = "b3",
  alive = "b5",
  age = "b7",
  age.truncate = 24,
  date.interview = "v008",
  month.cut = c(1, 12, 24, 36, 48, 60),
  year.cut = seq(1980, 2020, by = 5),
  min.last.period = 0,
  cmc.adjust = 0,
  compact = FALSE,
  compact.by = c("v001", "v024", "v025", "v005")
)

Arguments

Value

This function returns a new data frame where each row indicate a person-month, with the additional variables specified in the function argument.

Author(s)

Zehang Richard Li, Bryan Martin, Laina Mercer

References

Li, Z., Hsiao, Y., Godwin, J., Martin, B. D., Wakefield, J., Clark, S. J., & with support from the United Nations Inter-agency Group for Child Mortality Estimation and its technical advisory group. (2019). Changes in the spatial distribution of the under-five mortality rate: Small-area analysis of 122 DHS surveys in 262 subregions of 35 countries in Africa. PloS one, 14(1), e0210645.

Mercer, L. D., Wakefield, J., Pantazis, A., Lutambi, A. M., Masanja, H., & Clark, S. (2015). Space-time smoothing of complex survey data: small area estimation for child mortality. The annals of applied statistics, 9(4), 1889.

Examples

## Not run: 
my_fp <- "/myExampleFilepath/surveyData.DTA"
DemoData <- getBirths(filepath = my_fp, surveyyear = 2015) 

## End(Not run)

Description

Aggregate person-month data into counts and totals by groups.

Usage

getCounts(data, variables, by, ignore = NULL, addtotal = TRUE, drop = TRUE)

Arguments

Value

data.frame of the ggregated counts.

Author(s)

Zehang Richard Li

Examples

# a toy dataset with 4 time periods but one missing in data
timelist <- factor(1:4)
data = data.frame(died = c(0,0,0,1,1,0,0), 
					area = c(rep(c("A", "B"), 3), "A"), 
					time = timelist[c(1,1,2,3,3,3,3)])
data
# without ignore argument, all levels will be imputed
getCounts(data, variables = "died", by = c("area", "time"))

# ignoring time = 4, the ignored level will not be imputed (but still in the output)
getCounts(data, variables = "died", by = c("area", "time"), 
			ignore = list("time"=c(4)) )

Description

Extract posterior summaries of random effects

Usage

getDiag(
  inla_mod,
  field = c("space", "time", "spacetime")[1],
  CI = 0.95,
  draws = NULL,
  nsim = 1000,
  ...
)

Arguments

Value

List of diagnostic plots

Author(s)

Zehang Richard Li

Examples

## Not run: 
  data(DemoMap)
  years <- levels(DemoData[[1]]$time)
  
  # obtain direct estimates
  data <- getDirectList(births = DemoData, 
  years = years,  
  regionVar = "region", timeVar = "time", 
  clusterVar = "~clustid+id", 
  ageVar = "age", weightsVar = "weights", 
  geo.recode = NULL)
  # obtain direct estimates
  data_multi <- getDirectList(births = DemoData, years = years,
    regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
    ageVar = "age", weightsVar = "weights", geo.recode = NULL)
  data <- aggregateSurvey(data_multi)
  
  #  national model
  years.all <- c(years, "15-19")
  fit1 <- smoothDirect(data = data, geo = DemoMap$geo, Amat = DemoMap$Amat, 
    year_label = years.all, year_range = c(1985, 2019), 
    rw = 2, is.yearly=FALSE, m = 5)
random.time <- getDiag(fit1, field = "time")
  random.space <- getDiag(fit1, field = "space")
  random.spacetime <- getDiag(fit1, field = "spacetime")

## End(Not run)

Description

Obtain the Horvitz-Thompson direct estimates and standard errors using delta method for a single survey.

Usage

getDirect(
  births,
  years,
  regionVar = "region",
  timeVar = "time",
  clusterVar = "~v001+v002",
  ageVar = "age",
  weightsVar = "v005",
  Ntrials = NULL,
  geo.recode = NULL,
  national.only = FALSE,
  CI = 0.95
)

Arguments

Value

a matrix of period-region summary of the Horvitz-Thompson direct estimates by region and time period specified in the argument, the standard errors using delta method for a single survey, the 95% confidence interval, and the logit of the estimates.

Author(s)

Zehang Richard Li, Bryan Martin, Laina Mercer

References

Li, Z., Hsiao, Y., Godwin, J., Martin, B. D., Wakefield, J., Clark, S. J., & with support from the United Nations Inter-agency Group for Child Mortality Estimation and its technical advisory group. (2019). Changes in the spatial distribution of the under-five mortality rate: Small-area analysis of 122 DHS surveys in 262 subregions of 35 countries in Africa. PloS one, 14(1), e0210645.

Mercer, L. D., Wakefield, J., Pantazis, A., Lutambi, A. M., Masanja, H., & Clark, S. (2015). Space-time smoothing of complex survey data: small area estimation for child mortality. The annals of applied statistics, 9(4), 1889.

See Also

getDirectList

Examples

## Not run: 
data(DemoData)
years <- c("85-89", "90-94", "95-99", "00-04", "05-09", "10-14")
mean <- getDirect(births = DemoData[[1]],  years = years, 
regionVar = "region", timeVar = "time", clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", geo.recode = NULL)

## End(Not run)

Description

Obtain the Horvitz-Thompson direct estimates and standard errors using delta method for multiple surveys.

Usage

getDirectList(
  births,
  years,
  regionVar = "region",
  timeVar = "time",
  clusterVar = "~v001+v002",
  ageVar = "age",
  weightsVar = "v005",
  Ntrials = NULL,
  geo.recode = NULL,
  national.only = FALSE
)

Arguments

Value

This is the extension to the getDirect function that returns estimates from multiple surveys. Additional columns in the output (survey and surveyYears) specify the estimates from different surveys.

Author(s)

Zehang Richard Li, Bryan Martin, Laina Mercer

References

Li, Z., Hsiao, Y., Godwin, J., Martin, B. D., Wakefield, J., Clark, S. J., & with support from the United Nations Inter-agency Group for Child Mortality Estimation and its technical advisory group. (2019). Changes in the spatial distribution of the under-five mortality rate: Small-area analysis of 122 DHS surveys in 262 subregions of 35 countries in Africa. PloS one, 14(1), e0210645.

Mercer, L. D., Wakefield, J., Pantazis, A., Lutambi, A. M., Masanja, H., & Clark, S. (2015). Space-time smoothing of complex survey data: small area estimation for child mortality. The annals of applied statistics, 9(4), 1889.

See Also

getDirect

Examples

## Not run: 
data(DemoData)
years <- c("85-89", "90-94", "95-99", "00-04", "05-09", "10-14")
mean <- getDirectList(births = DemoData, years = years, 
regionVar = "region", timeVar = "time", clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", geo.recode = NULL)

## End(Not run)

Description

Extract smoothed estimates.

Usage

getSmoothed(
  inla_mod,
  nsim = 1000,
  weight.strata = NULL,
  weight.frame = NULL,
  verbose = FALSE,
  mc = 0,
  include_time_unstruct = FALSE,
  CI = 0.95,
  draws = NULL,
  save.draws = FALSE,
  include_subnational = TRUE,
  only.age = NULL,
  joint = FALSE,
  ...
)

Arguments

Value

A data frame or a list of data frames of S3 class SUMMERproj, which contains the smoothed estimates.

Author(s)

Zehang Richard Li

See Also

plot.SUMMERproj

Examples

## Not run: 
years <- levels(DemoData[[1]]$time)

# obtain direct estimates
data <- getDirectList(births = DemoData, 
years = years,  
regionVar = "region", timeVar = "time", 
clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", 
geo.recode = NULL)
# obtain direct estimates
data_multi <- getDirectList(births = DemoData, years = years,
  regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
  ageVar = "age", weightsVar = "weights", geo.recode = NULL)
data <- aggregateSurvey(data_multi)

##
## The following example shows extracting estimates from 
##   fitted smoothDirect() model
##
#  national model
years.all <- c(years, "15-19")
fit1 <- smoothDirect(data = data, Amat = NULL, 
  year_label = years.all, year_range = c(1985, 2019), 
  rw = 2, is.yearly=FALSE, m = 5)
out1 <- getSmoothed(fit1)
plot(out1, is.subnational=FALSE)

#  subnational model
fit2 <- smoothDirect(data = data, Amat = mat, 
  year_label = years.all, year_range = c(1985, 2019), 
  rw = 2, is.yearly=TRUE, m = 5, type.st = 4)
out2 <- getSmoothed(fit2)
plot(out2, is.yearly=TRUE, is.subnational=TRUE)

## The following examples compare the marginal posterior draws
##    with joint posterior draws. The latter gives draws of 
##    all estimates in addition to marginal summaries.
## The plots should fall closely along the 45 degree line 

#  national period model
years.all <- c(years, "15-19")
fit0 <- smoothDirect(data = data, Amat = NULL, 
year_label = years, year_range = c(1985, 2019), 
time.model = 'rw2', m = 1, control.compute = list(config =TRUE))
out0 <- getSmoothed(fit0)
out0a <- getSmoothed(fit0, joint = TRUE, nsim = 1e5)
par(mfrow = c(1, 3))
plot(out0$median, out0a$overall$median)
abline(c(0, 1))
plot(out0$upper, out0a$overall$upper)
abline(c(0, 1))
plot(out0$lower, out0a$overall$lower)
abline(c(0, 1))

#  national yearly model
years.all <- c(years, "15-19")
fit1 <- smoothDirect(data = data, Amat = NULL, 
year_label = years.all, year_range = c(1985, 2019), 
time.model = 'rw2', m = 5, control.compute = list(config =TRUE))
out1 <- getSmoothed(fit1)
out1a <- getSmoothed(fit1, joint = TRUE, nsim = 1e5, save.draws = TRUE)
par(mfrow = c(1, 3))
plot(out1$median, out1a$overall$median)
abline(c(0, 1))
plot(out1$upper, out1a$overall$upper)
abline(c(0, 1))
plot(out1$lower, out1a$overall$lower)
abline(c(0, 1))

#  subnational model
fit2 <- smoothDirect(data = data, Amat = DemoMap$Amat, 
year_label = years.all, year_range = c(1985, 2019), 
time.model = 'rw2',is.yearly=TRUE, m = 5, type.st = 4, 
control.compute = list(config =TRUE))
out2 <- getSmoothed(fit2)
out2a <- getSmoothed(fit2, joint = TRUE, nsim = 1e5, save.draws = TRUE)
par(mfrow = c(1, 3))
plot(out2$median, out2a$overall$median)
abline(c(0, 1))
plot(out2$upper, out2a$overall$upper)
abline(c(0, 1))
plot(out2$lower, out2a$overall$lower)
abline(c(0, 1))

#  subnational space-only model combining all periods
fit3 <- smoothDirect(data = data, 
       time.model = NULL, Amat = DemoMap$Amat, 
   control.compute = list(config =TRUE))
out3 <- getSmoothed(fit3)
out3a <- getSmoothed(fit3, joint = TRUE, nsim = 1e5, save.draws = TRUE)
par(mfrow = c(1, 3))
plot(out3$median, out3a$overall$median)
abline(c(0, 1))
plot(out3$upper, out3a$overall$upper)
abline(c(0, 1))
plot(out3$lower#' , out3a$overall$lower)
abline(c(0, 1))

## End(Not run)

Description

This function visualizes the map with different variables. The input data frame can be either the long or wide format.

Usage

hatchPlot(
  data,
  variables,
  values = NULL,
  labels = NULL,
  geo,
  by.data,
  by.geo,
  is.long = FALSE,
  lower,
  upper,
  lim = NULL,
  lim.CI = NULL,
  breaks.CI = NULL,
  ncol = 4,
  hatch = NULL,
  border = NULL,
  size = 1,
  legend.label = NULL,
  per1000 = FALSE,
  direction = 1,
  ...
)

Arguments

Author(s)

Zehang Richard Li, Katie Wilson

Examples

## Not run: 
years <- levels(DemoData[[1]]$time)

# obtain direct estimates
data <- getDirectList(births = DemoData, 
years = years,  
regionVar = "region", timeVar = "time", 
clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", 
geo.recode = NULL)
# obtain direct estimates
data_multi <- getDirectList(births = DemoData, years = years,
  regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
  ageVar = "age", weightsVar = "weights", geo.recode = NULL)
data <- aggregateSurvey(data_multi)

fit2 <- smoothDirect(data = data, geo = geo, Amat = mat, 
  year_label = years.all, year_range = c(1985, 2019), 
  rw = 2, is.yearly=TRUE, m = 5, type.st = 4)
out2 <- getSmoothed(fit2)

plot(out2, is.yearly=TRUE, is.subnational=TRUE)

hatchPlot(data = subset(out2, is.yearly==FALSE), geo = geo,
variables=c("years"), values = c("median"), 
by.data = "region", by.geo = "REGNAME", 
lower = "lower", upper = "upper", is.long=TRUE)


## End(Not run)

Description

New random IID models for m-year to period random effects

Usage

iid.new(
  cmd = c("graph", "Q", "mu", "initial", "log.norm.const", "log.prior", "quit"),
  theta = NULL
)

Arguments

Description

New random IID models for m-year to period random effects

Usage

iid.new.pc(
  cmd = c("graph", "Q", "mu", "initial", "log.norm.const", "log.prior", "quit"),
  theta = NULL
)

Arguments

Description

The list contains several data frames.

Usage

data(KenData)

Format

An object of class list of length 4.

Details

• HIV2014, a data frame with three columns: years (in five year periods), region (8 Admin-1 region groups), and the estimated bias of the reported U5MR due to HIV for each 5 year period from 1990-1994 to 2010-2014. The bias is represented as the ratio of the reported U5MR to the true U5MR.

• HIV2014.yearly, a data frame with three columns: years (in one year interval), region (8 Admin-1 region groups), and the estimated bias of the reported U5MR due to HIV for each year from 1980 to 2014. The bias is represented as the ratio of the reported U5MR to the true U5MR.

• IGME2019. Yearly Estimates of national under-5 child mortality in Kenya from the 2019 UN-IGME estimates.

• UrbanProp. Proportion of urban population by county and total population by county. Source: 2009 Kenya Population and Housing Census, and Table A2 of Kenya 2014 DHS report.

References

Neff Walker, Kenneth Hill, and Fengmin Zhao (2012) Child mortality estimation: methods used to adjust for bias due to aids in estimating trends in under-five mortality.,
PLoS Medicine, 9(8):e1001298.

Description

Datasets related to the 2009 census frame for Kenya based on the 2009 Kenya Population and Housing Census. General population totals are estimated for 2014. Based on 2014 population density estimates interpolated with exponential growth rate between 2010 and 2015 from WorldPop data.

Usage

data(kenyaPopulationData)

easpaKenyaNeonatal

poppaKenya

poppsubKenya

Format

A number of data.frames with information about the 2009 Kenya Population and Housing Census and the population in Kenya at the time of the 2014 Demographic Health Survey. Some of the data.frames have been adjusted to contain information about neonatals born from 2010-2014 rather than general population in 2014. The dataset names are: easpaKenya, easpaKenyaNeonatal, poppaKenya, and poppsubKenya.

An object of class data.frame with 47 rows and 12 columns.

An object of class data.frame with 47 rows and 6 columns.

An object of class data.frame with 300 rows and 7 columns.

Source

<http://dhsprogram.com/pubs/pdf/FR308/FR308.pdf>

References

Kenya National Bureau of Statistics, Ministry of Health/Kenya, National AIDS Control Council/Kenya, Kenya Medical Research Institute, and National Council For Population And Development/Kenya, 2015. Kenya Demographic and Health Survey 2014. Rockville, Maryland, USA. URL: http://dhsprogram.com/pubs/pdf/FR308/FR308.pdf.

Stevens, F.R., Gaughan, A.E., Linard, C., Tatem, A.J., 2015. Disaggregat- ing census data for population mapping using random forests with remotely-sensed and ancillary data. PloS One 10, e0107042.

Tatem, A.J., 2017. WorldPop, open data for spatial demography. Scientific Data 4.

Description

Shapefiles are King County in the Washington States.

Usage

KingCounty

Format

An object of class SpatialPolygonsDataFrame with 48 rows and 9 columns.

Description

Logit transformation

Usage

logit(x)

Arguments

Value

logit of x

Examples

x <- .5
logit(x)

Description

Adapted from logitnorm package. Calculates the mean of a distribution whose logit is Gaussian. Each row of muSigmaMat is a mean and standard deviation on the logit scale.

Usage

logitNormMean(muSigmaMat, logisticApprox = FALSE, ...)

Arguments

Details

If $\mbox{logit}(Y) \sim N(\mu, \sigma^2)$, This function calculates $E[Y]$ via either numerical integration or by assuming that Y follows a logistic distribution. Under this approximation, setting $k = 16 \sqrt(3) / (15 \pi)$, we approximate the expectation as:

$E[Y] = expit(\mu / \sqrt(1 + k^2 \sigma^2))$

. The above logistic approximation speeds up the computation, but also sacrifices some accuracy.

Value

A vector of expectations of the specified random variables

Author(s)

John Paige

Examples

mus = c(-5, 0, 5)
sigmas = rep(1, 3)
logitNormMean(cbind(mus, sigmas))
logitNormMean(cbind(mus, sigmas), TRUE)

Description

Functions for generating pixellated population information and population frames at the 'area' and 'subarea' levels. The 'area' and 'subarea' levels can be thought of as big regions and little regions, where areas can be partitioned into unique sets of subareas. For example, Admin-1 and Admin-2 areas might be areas and subareas respectively. The population totals are either tabulated at the area x urban/rural level, the subarea x urban/rural level, or at the pixel level of a specified resolution. Totals are calculated using population density information, shapefiles, and, possibly, preexisting population frames at different areal levels. Note that area names should each be unique, and similarly for subarea names.

Usage

makePopIntegrationTab(
  kmRes = 5,
  pop,
  domainMapDat,
  eastLim,
  northLim,
  mapProjection,
  areaMapDat,
  subareaMapDat,
  areaNameVar = "NAME_1",
  subareaNameVar = "NAME_2",
  poppa = NULL,
  poppsub = NULL,
  stratifyByUrban = TRUE,
  areapa = NULL,
  areapsub = NULL,
  customSubsetPolygons = NULL,
  areaPolygonSubsetI = NULL,
  subareaPolygonSubsetI = NULL,
  mean.neighbor = 50,
  delta = 0.1,
  returnPoppTables = FALSE,
  setNAsToZero = TRUE,
  fixZeroPopDensitySubareas = FALSE,
  extractMethod = "bilinear"
)

getPoppsub(
  kmRes = 1,
  pop,
  domainMapDat,
  eastLim,
  northLim,
  mapProjection,
  poppa,
  areapa = NULL,
  areapsub,
  subareaMapDat,
  subareaNameVar = "NAME_2",
  stratifyByUrban = TRUE,
  areaMapDat = NULL,
  areaNameVar = "NAME_1",
  areaPolygonSubsetI = NULL,
  subareaPolygonSubsetI = NULL,
  customSubsetPolygons = NULL,
  mean.neighbor = 50,
  delta = 0.1,
  setNAsToZero = TRUE,
  fixZeroPopDensitySubareas = FALSE
)

adjustPopMat(
  popMat,
  poppaTarget = NULL,
  adjustBy = c("area", "subarea"),
  stratifyByUrban = TRUE
)

Arguments

Functions

• makePopIntegrationTab(): Generate pixellated 'grid' of coordinates (both longitude/latitude and east/north) over spatial domain of the given resolution with associated population totals, areas, subareas, and urban/rural levels. For very small areas that might not otherwise have a grid point in them, a custom integration point is added at their centroid. Sets urbanicity classifications by thresholding input population density raster using area and subarea population tables, and generates area and subarea population tables from population density information if not already given. Can be used for integrating predictions from the given coordinates to area and subarea levels using population weights.

• getPoppsub(): Generate table of estimates of population totals per subarea x urban/rural combination based on population density raster at 'kmres' resolution "grid", including custom integration points for any subarea too small to include grid points at their centroids.

• adjustPopMat(): Adjust population densities in grid based on a population frame.

Author(s)

John Paige

See Also

setThresholdsByRegion, poppRegionFromPopMat, simPopSPDE, simPopCustom

Examples

## Not run: 

library(sp)
library(sf)
# download Kenya GADM shapefiles from SUMMERdata github repository
githubURL <- paste0("https://github.com/paigejo/SUMMERdata/blob/main/data/", 
                    "kenyaMaps.rda?raw=true")
tempDirectory = "~/"
mapsFilename = paste0(tempDirectory, "/kenyaMaps.rda")
if(!file.exists(mapsFilename)) {
  download.file(githubURL,mapsFilename)
}

# load it in
out = load(mapsFilename)
out
adm1@data$NAME_1 = as.character(adm1@data$NAME_1)
adm1@data$NAME_1[adm1@data$NAME_1 == "Trans Nzoia"] = "Trans-Nzoia"
adm1@data$NAME_1[adm1@data$NAME_1 == "Elgeyo-Marakwet"] = "Elgeyo Marakwet"
adm2@data$NAME_1 = as.character(adm2@data$NAME_1)
adm2@data$NAME_1[adm2@data$NAME_1 == "Trans Nzoia"] = "Trans-Nzoia"
adm2@data$NAME_1[adm2@data$NAME_1 == "Elgeyo-Marakwet"] = "Elgeyo Marakwet"

# some Admin-2 areas have the same name
adm2@data$NAME_2 = as.character(adm2@data$NAME_2)
adm2@data$NAME_2[(adm2@data$NAME_1 == "Bungoma") & 
                   (adm2@data$NAME_2 == "Lugari")] = "Lugari, Bungoma"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Kakamega") & 
                   (adm2@data$NAME_2 == "Lugari")] = "Lugari, Kakamega"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Meru") & 
                   (adm2@data$NAME_2 == "Igembe South")] = "Igembe South, Meru"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Tharaka-Nithi") & 
                   (adm2@data$NAME_2 == "Igembe South")] = "Igembe South, Tharaka-Nithi"

# The spatial area of unknown 8 is so small, it causes problems unless its removed or 
# unioned with another subarea. Union it with neighboring Kakeguria:
newadm2 = adm2
unknown8I = which(newadm2$NAME_2 == "unknown 8")
newadm2$NAME_2[newadm2$NAME_2 %in% c("unknown 8", "Kapenguria")] <- 
  "Kapenguria + unknown 8"
admin2.IDs <- newadm2$NAME_2

newadm2@data = cbind(newadm2@data, NAME_2OLD = newadm2@data$NAME_2)
newadm2@data$NAME_2OLD = newadm2@data$NAME_2
newadm2@data$NAME_2 = admin2.IDs
newadm2$NAME_2 = admin2.IDs
temp <- terra::aggregate(as(newadm2, "SpatVector"), by="NAME_2")

temp <- sf::st_as_sf(temp)
temp <- sf::as_Spatial(temp)

tempData = newadm2@data[-unknown8I,]
tempData = tempData[order(tempData$NAME_2),]
newadm2 <- SpatialPolygonsDataFrame(temp, tempData, match.ID = F)
adm2 = newadm2

# download 2014 Kenya population density TIF file

githubURL <- paste0("https://github.com/paigejo/SUMMERdata/blob/main/data/", 
                    "Kenya2014Pop/worldpop_total_1y_2014_00_00.tif?raw=true")
popTIFFilename = paste0(tempDirectory, "/worldpop_total_1y_2014_00_00.tif")
if(!file.exists(popTIFFilename)) {
  download.file(githubURL,popTIFFilename)
}

# load it in
pop = terra::rast(popTIFFilename)

eastLim = c(-110.6405, 832.4544)
northLim = c(-555.1739, 608.7130)

## Construct poppsubKenya, a table of urban/rural general population totals 
## in each subarea. Technically, this is not necessary since we can load in 
## poppsubKenya via data(kenyaPopulationData). First, we will need to calculate 
## the areas in km^2 of the areas and subareas


# use Lambert equal area projection of areas (Admin-1) and subareas (Admin-2)
midLon = mean(adm1@bbox[1,])
midLat = mean(adm1@bbox[2,])
p4s = paste0("+proj=laea +x_0=0 +y_0=0 +lon_0=", midLon, 
             " +lat_0=", midLat, " +units=km")

adm1_sf = st_as_sf(adm1)
adm1proj_sf = st_transform(adm1_sf, p4s)
adm1proj = as(adm1proj_sf, "Spatial")

adm2_sf = st_as_sf(adm2)
adm2proj_sf = st_transform(adm2_sf, p4s)
adm2proj = as(adm2proj_sf, "Spatial")

# now calculate spatial area in km^2
admin1Areas = as.numeric(st_area(adm1proj_sf))
admin2Areas = as.numeric(st_area(adm2proj_sf))

areapaKenya = data.frame(area=adm1proj@data$NAME_1, spatialArea=admin1Areas)
areapsubKenya = data.frame(area=adm2proj@data$NAME_1, subarea=adm2proj@data$NAME_2, 
                           spatialArea=admin2Areas)

# Calculate general population totals at the subarea (Admin-2) x urban/rural 
# level and using 1km resolution population grid. Assign urbanicity by 
# thresholding population density based on estimated proportion population 
# urban/rural, making sure total area (Admin-1) urban/rural populations in 
# each area matches poppaKenya.
require(fields)
# NOTE: the following function will typically take ~15-20 minutes. Can speed up 
#       by setting kmRes to be higher, but we recommend fine resolution for 
#       this step, since it only needs to be done once. Instead of running this, 
#       you can simply run data(kenyaPopulationData)
system.time(poppsubKenya <- getPoppsub(
  kmRes=1, pop=pop, domainMapDat=adm0,
  eastLim=eastLim, northLim=northLim, mapProjection=projKenya,
  poppa = poppaKenya, areapa=areapaKenya, areapsub=areapsubKenya, 
  areaMapDat=adm1, subareaMapDat=adm2, 
  areaNameVar = "NAME_1", subareaNameVar="NAME_2"))

# Now generate a general population integration table at 5km resolution, 
# based on subarea (Admin-2) x urban/rural population totals. This takes 
# ~1 minute
system.time(popMatKenya <- makePopIntegrationTab(
  kmRes=5, pop=pop, domainMapDat=adm0,
  eastLim=eastLim, northLim=northLim, mapProjection=projKenya,
  poppa = poppaKenya, poppsub=poppsubKenya, 
  areaMapDat = adm1, subareaMapDat = adm2,
  areaNameVar = "NAME_1", subareaNameVar="NAME_2"))

## Adjust popMat to be target (neonatal) rather than general population density. First
## creat the target population frame
## (these numbers are based on IPUMS microcensus data)
mothersPerHouseholdUrb = 0.3497151
childrenPerMotherUrb = 1.295917
mothersPerHouseholdRur = 0.4787696
childrenPerMotherRur = 1.455222
targetPopPerStratumUrban = easpaKenya$HHUrb * mothersPerHouseholdUrb * childrenPerMotherUrb
targetPopPerStratumRural = easpaKenya$HHRur * mothersPerHouseholdRur * childrenPerMotherRur
easpaKenyaNeonatal = easpaKenya
easpaKenyaNeonatal$popUrb = targetPopPerStratumUrban
easpaKenyaNeonatal$popRur = targetPopPerStratumRural
easpaKenyaNeonatal$popTotal = easpaKenyaNeonatal$popUrb + easpaKenyaNeonatal$popRur
easpaKenyaNeonatal$pctUrb = 100 * easpaKenyaNeonatal$popUrb / easpaKenyaNeonatal$popTotal
easpaKenyaNeonatal$pctTotal = 
  100 * easpaKenyaNeonatal$popTotal / sum(easpaKenyaNeonatal$popTotal)

# Generate the target population density by scaling the current population density grid 
# at the Admin1 x urban/rural level
popMatKenyaNeonatal = adjustPopMat(popMatKenya, easpaKenyaNeonatal)

# Generate neonatal population table from the neonatal population integration matrix.
# This is technically not necessary for population simulation purposes, but is here 
# for illustrative purposes
poppsubKenyaNeonatal = poppRegionFromPopMat(popMatKenyaNeonatal, popMatKenyaNeonatal$subarea)
poppsubKenyaNeonatal = cbind(subarea=poppsubKenyaNeonatal$region, 
                             area=adm2@data$NAME_1[match(poppsubKenyaNeonatal$region, 
                               adm2@data$NAME_2)], 
                             poppsubKenyaNeonatal[,-1])
print(head(poppsubKenyaNeonatal))

## End(Not run)

Description

The list contains several data frames.

Usage

data(MalawiData)

Format

An object of class list of length 4.

Details

• HIV, a data frame with three columns: years (in five year periods), survey, and the estimated bias of the reported U5MR due to HIV for each 5 year period. The bias is represented as the ratio of the reported U5MR to the true U5MR.

• HIV.yearly, a data frame with three columns: years (in one year interval), survey, and the estimated bias of the reported U5MR due to HIV for each year. The bias is represented as the ratio of the reported U5MR to the true U5MR.

• IGME2019. Yearly Estimates of national under-5 child mortality in Malawi from the 2019 UN-IGME estimates.

• IGME2019.nmr. Yearly Estimates of national neonatal mortality in Malawi from the 2019 UN-IGME estimates.

References

Neff Walker, Kenneth Hill, and Fengmin Zhao (2012) Child mortality estimation: methods used to adjust for bias due to aids in estimating trends in under-five mortality.,
PLoS Medicine, 9(8):e1001298.

Description

SpatialPolygonsDataFrame objects that reflect the Admin 2 regions in Malawi, including the Likoma island. The Admin 2 region names are in the ADM2_EN field.

Usage

data(MalawiMap)

Format

An object of class SpatialPolygonsDataFrame with 28 rows and 14 columns.

Description

Mapping estimates for svysae object

Usage

mapEstimates(x, geo.data, variable, viridis.option = "viridis")

Arguments

Value

ggplot containing map of small area posterior summary statistics

Examples

## Not run: 
data(DemoData2)
data(DemoMap2)
library(survey)
des0 <- svydesign(ids = ~clustid+id, strata = ~strata,
                  weights = ~weights, data = DemoData2, nest = TRUE)
Xmat <- aggregate(age~region, data = DemoData2, FUN = mean)
geo.data <- sf::st_as_sf(DemoMap2$geo)
geo.data$domain <- geo.data$REGNAME
cts.res <- smoothArea(tobacco.use ~ 1,
                      domain = ~region,
                      design = des0,
                      adj.mat = DemoMap2$Amat, 
                      pc.u = 1,
                      pc.alpha = 0.01,
                      pc.u.phi = 0.5,
                      pc.alpha.phi = 2/3,
                      return.samples = TRUE)
mapEstimates(cts.res, geo.data = geo.data, variable = "median")
mapEstimates(cts.res, geo.data = geo.data, variable = "var")

## End(Not run)

Description

This function visualizes the map with different variables. The input data frame can be either the long or wide format.

Usage

mapPlot(
  data = NULL,
  variables,
  values = NULL,
  labels = NULL,
  geo,
  by.data,
  by.geo,
  is.long = FALSE,
  size = 0.5,
  removetab = FALSE,
  border = "gray20",
  ncol = NULL,
  ylim = NULL,
  legend.label = NULL,
  per1000 = FALSE,
  clean = TRUE,
  size.label = 2,
  add.adj = FALSE,
  color.adj = "red",
  alpha.adj = 0.85,
  direction = 1,
  cut = NULL
)

Arguments

Author(s)

Zehang Richard Li

Examples

## Not run: 
data(DemoMap)
# Plotting data in the long format
dat <- data.frame(region = rep(c("central",  "eastern", "northern", "western"), 3),
year = rep(c(1980, 1990, 2000), each = 4),
values = stats::rnorm(12))
utils::head(dat)
mapPlot(dat, variables = "year", values = "values",
by.data = "region", geo = DemoMap$geo,
by.geo = "NAME_final", is.long = TRUE)
dat <- data.frame(region = c("central",  "eastern", "northern", "western"),
Year1 = stats::rnorm(4), Year2 = stats::rnorm(4),
Year3 = stats::rnorm(4))
utils::head(dat)
mapPlot(dat, variables = c("Year1", "Year2", "Year3"),
 labels = c(1980, 1990, 2000),
by.data = "region", geo = DemoMap$geo,
by.geo = "NAME_final", is.long = FALSE)


## End(Not run)

Description

Map GPS points to polygon regions

Usage

mapPoints(data, geo, long, lat, names)

Arguments

Value

Spatial djacency matrix.

Author(s)

Zehang Richard Li

Examples

data(DemoMap) 
dat <- data.frame(ID = c(1,2,3), lon = c(32.2, 33.7, 33), lat = c(0.1, 0.9, 2.8))
dat2 <- mapPoints(dat, DemoMap$geo, long = "lon", lat = "lat", names = "REGNAME")
dat2

Description

Plot projection output.

Usage

## S3 method for class 'SUMMERproj'
plot(
  x,
  year_label = c("85-89", "90-94", "95-99", "00-04", "05-09", "10-14", "15-19"),
  year_med = c(1987, 1992, 1997, 2002, 2007, 2012, 2017),
  is.subnational = TRUE,
  proj_year = 2015,
  data.add = NULL,
  option.add = list(point = NULL, lower = NULL, upper = NULL, by = NULL),
  color.add = "black",
  label.add = NULL,
  dodge.width = 1,
  plot.CI = NULL,
  per1000 = FALSE,
  color.CI = NULL,
  alpha.CI = 0.5,
  ...
)

Arguments

Author(s)

Zehang Richard Li

See Also

getSmoothed

Examples

## Not run: 
years <- levels(DemoData[[1]]$time)

# obtain direct estimates
data <- getDirectList(births = DemoData, 
years = years,  
regionVar = "region", timeVar = "time", 
clusterVar = "~clustid+id", 
ageVar = "age", weightsVar = "weights", 
geo.recode = NULL)
# obtain direct estimates
data_multi <- getDirectList(births = DemoData, years = years,
  regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
  ageVar = "age", weightsVar = "weights", geo.recode = NULL)
data <- aggregateSurvey(data_multi)

#  national model
years.all <- c(years, "15-19")
fit1 <- smoothDirect(data = data, geo = NULL, Amat = NULL, 
  year_label = years.all, year_range = c(1985, 2019), 
  rw = 2, is.yearly=FALSE, m = 5)
out1 <- getSmoothed(fit1)
plot(out1, is.subnational=FALSE)

#  subnational model
fit2 <- smoothDirect(data = data, geo = geo, Amat = mat, 
  year_label = years.all, year_range = c(1985, 2019), 
  rw = 2, is.yearly=TRUE, m = 5, type.st = 4)
out2 <- getSmoothed(fit2)
plot(out2, is.yearly=TRUE, is.subnational=TRUE)



## End(Not run)

Generate a population frame of a similar format to poppa argument of simPopCustom with a custom set of regions

Description

Generate a population frame of a similar format to poppa argument of simPopCustom with a custom set of regions

Usage

poppRegionFromPopMat(popMat, regions)

Arguments

Details

Urbanicity thresholds are set based on that region's percent population urban. Intended as a helper function of getPoppsub, but can also be used for custom sets of regions (i.e. more than just 2 areal levels: area and subarea).

Value

A table of population totals by region

Author(s)

John Paige

See Also

getPoppsub

Examples

## Not run: 
data(kenyaPopulationData)

#' # download Kenya GADM shapefiles from SUMMERdata github repository
githubURL <- "https://github.com/paigejo/SUMMERdata/blob/main/data/kenyaMaps.rda?raw=true"
tempDirectory = "~/"
mapsFilename = paste0(tempDirectory, "/kenyaMaps.rda")
if(!file.exists(mapsFilename)) {
  download.file(githubURL,mapsFilename)
}

# load it in
out = load(mapsFilename)
out
adm1@data$NAME_1 = as.character(adm1@data$NAME_1)
adm1@data$NAME_1[adm1@data$NAME_1 == "Trans Nzoia"] = "Trans-Nzoia"
adm1@data$NAME_1[adm1@data$NAME_1 == "Elgeyo-Marakwet"] = "Elgeyo Marakwet"
adm2@data$NAME_1 = as.character(adm2@data$NAME_1)
adm2@data$NAME_1[adm2@data$NAME_1 == "Trans Nzoia"] = "Trans-Nzoia"
adm2@data$NAME_1[adm2@data$NAME_1 == "Elgeyo-Marakwet"] = "Elgeyo Marakwet"

# some Admin-2 areas have the same name
adm2@data$NAME_2 = as.character(adm2@data$NAME_2)
adm2@data$NAME_2[(adm2@data$NAME_1 == "Bungoma") & 
  (adm2@data$NAME_2 == "Lugari")] = "Lugari, Bungoma"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Kakamega") & 
  (adm2@data$NAME_2 == "Lugari")] = "Lugari, Kakamega"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Meru") & 
  (adm2@data$NAME_2 == "Igembe South")] = "Igembe South, Meru"
adm2@data$NAME_2[(adm2@data$NAME_1 == "Tharaka-Nithi") & 
  (adm2@data$NAME_2 == "Igembe South")] = "Igembe South, Tharaka-Nithi"

# The spatial area of unknown 8 is so small, it causes problems unless 
# its removed or unioned with another subarea. Union it with neighboring 
# Kakeguria:
newadm2 = adm2
unknown8I = which(newadm2$NAME_2 == "unknown 8")
newadm2$NAME_2[newadm2$NAME_2 %in% c("unknown 8", "Kapenguria")] <- "Kapenguria + unknown 8"
admin2.IDs <- newadm2$NAME_2

newadm2@data = cbind(newadm2@data, NAME_2OLD = newadm2@data$NAME_2)
newadm2@data$NAME_2OLD = newadm2@data$NAME_2
newadm2@data$NAME_2 = admin2.IDs
newadm2$NAME_2 = admin2.IDs
temp <- terra::aggregate(as(newadm2, "SpatVector"), by="NAME_2")

library(sf)
temp <- sf::st_as_sf(temp)
temp <- sf::as_Spatial(temp)

tempData = newadm2@data[-unknown8I,]
tempData = tempData[order(tempData$NAME_2),]
newadm2 <- SpatialPolygonsDataFrame(temp, tempData, match.ID = F)
adm2 = newadm2

# download 2014 Kenya population density TIF file

githubURL <- paste0("https://github.com/paigejo/SUMMERdata/blob/main/data/", 
                    "Kenya2014Pop/worldpop_total_1y_2014_00_00.tif?raw=true")
popTIFFilename = paste0(tempDirectory, "/worldpop_total_1y_2014_00_00.tif")
if(!file.exists(popTIFFilename)) {
  download.file(githubURL,popTIFFilename)
}

# load it in
pop = terra::rast(popTIFFilename)

eastLim = c(-110.6405, 832.4544)
northLim = c(-555.1739, 608.7130)

require(fields)

# Now generate a general population integration table at 5km resolution, 
# based on subarea (Admin-2) x urban/rural population totals. This takes 
# ~1 minute
system.time(popMatKenya <- makePopIntegrationTab(
  kmRes=5, pop=pop, domainMapDat=adm0,
  eastLim=eastLim, northLim=northLim, mapProjection=projKenya,
  poppa = poppaKenya, poppsub=poppsubKenya, 
  areaMapDat = adm1, subareaMapDat = adm2,
  areaNameVar = "NAME_1", subareaNameVar="NAME_2"))
  
out = poppRegionFromPopMat(popMatKenya, popMatKenya$area)
out
poppaKenya

out = poppRegionFromPopMat(popMatKenya, popMatKenya$subarea)
out
poppsubKenya

popMatKenyaUnstratified = popMatKenya
popMatKenyaUnstratified$urban = NULL
out = poppRegionFromPopMat(popMatKenyaUnstratified, popMatKenyaUnstratified$area)
out
poppaKenya

## End(Not run)

Description

This function is the print method for class SUMMERmodel.

Usage

## S3 method for class 'SUMMERmodel'
print(x, ...)

Arguments

Author(s)

Zehang Li

See Also

summary.SUMMERmodel

Examples

## Not run: 
  library(SUMMER)
  library(dplyr)
  data(DemoData)

  # Smooth Direct Model
  years <- levels(DemoData[[1]]$time)
  # obtain direct estimates
  data_multi <- getDirectList(births = DemoData, years = years,
  regionVar = "region",  timeVar = "time", clusterVar = "~clustid+id",
  ageVar = "age", weightsVar = "weights", geo.recode = NULL)
  data <- aggregateSurvey(data_multi)
  
  years.all <- c(years, "15-19")
  fit <- smoothDirect(data = data, Amat = NULL, 
  year_label = years.all, year_range = c(1985, 2019), 
  time.model = 'rw2', is.yearly=FALSE, m = 5)
  fit

  # Cluster-level Model
  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"))
  fit

## End(Not run)

Description

This function is the print method for class SUMMERmodel.svy.

Usage

## S3 method for class 'SUMMERmodel.svy'
print(x, ...)

Arguments

Author(s)

Zehang Li

See Also

summary.SUMMERmodel.svy

Examples

## Not run: 
data(DemoData2)
data(DemoMap2)
fit0 <- smoothSurvey(data=DemoData2,  
Amat=DemoMap2$Amat, responseType="binary", 
responseVar="tobacco.use", strataVar="strata", 
weightVar="weights", regionVar="region", 
clusterVar = "~clustid+id", CI = 0.95)
fit0

## End(Not run)

Description

This function is the print method for class SUMMERprojlist.

Usage

## S3 method for class 'SUMMERprojlist'
print(x, ...)

Arguments

Author(s)

Zehang Li

Examples

## Not run: 
 library(SUMMER)
 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)

## End(Not run)

Description

Projection specifically chosen for Kenya. Project from lat/lon to northing/easting in kilometers. Uses epsg=21097 with km units. May not work on all systems due to differences in the behavior between different PROJ and GDAL versions.

Usage

projKenya(lon, lat = NULL, inverse = FALSE)

Arguments