Package {fastsae}

Type:	Package
Title:	Fast Implementation of Small Area Estimation Methods
Version:	0.1.0
Description:	Provides high-performance implementations of Small Area Estimation (SAE) methods leveraging C++ ('Rcpp', 'RcppArmadillo') and 'OpenMP' multi-threading. Supports standard area-level Fay-Herriot models (Fay and Herriot, 1979 <doi:10.1080/01621459.1979.10482505>), Spatial Fay-Herriot models (Pratesi and Salvati, 2008 <doi:10.1002/env.861>), Spatio-Temporal Fay-Herriot models (Marhuenda et al., 2013 <doi:10.1016/j.csda.2013.01.016>), and unit-level Battese-Harter-Fuller models (Battese et al., 1988 <doi:10.1080/01621459.1988.10478561>). Features include empirical best linear unbiased prediction (EBLUP), analytical and bootstrap Mean Squared Error (MSE) estimation, automatic handling of unsampled domains, and modern S3 diagnostic methods.
License:	GPL (≥ 3)
URL:	https://ridsonap.github.io/fastsae/, https://github.com/ridsonap/fastsae
BugReports:	https://github.com/ridsonap/fastsae/issues
Encoding:	UTF-8
LazyData:	true
Imports:	cli, ggplot2, lme4, methods, Rcpp, rlang, stats, utils
Suggests:	doParallel, dplyr, emdi, foreach, MASS, parallel, sae, knitr, rmarkdown, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
LinkingTo:	Rcpp, RcppArmadillo
Depends:	R (≥ 4.1.0)
Collate:	utils.R RcppExports.R data.R methods.R autoplot.R eblup_fh.R eblup_sfh.R eblup_stfh.R eblup_bhf.R
Config/roxygen2/version:	8.0.0
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2026-09-21 13:12:35 UTC; macbook
Author:	Ridson Al Farizal P [aut, cre, cph], Azka Ubaidillah [aut]
Maintainer:	Ridson Al Farizal P <alfrzlp@gmail.com>
Repository:	CRAN
Date/Publication:	2026-09-30 09:50:02 UTC

Autoplot Method for fastsae Objects

Description

Creates diagnostic and comparison plots for Small Area Estimation (SAE) models fitted with fastsae. This extends the generic autoplot from ggplot2.

Usage

## S3 method for class 'fastsae'
autoplot(object, type = c("comparison", "mse", "estimates", "scatter"), ...)

## S3 method for class 'list'
autoplot(object, type = c("comparison", "mse", "scatter"), ...)

Arguments

Value

A ggplot object.

Examples

library(fastsae)

# Single model plot
fit_fh <- eblup_fh(y ~ x1 + x2 + x3, data = mys, vardir = "vardir")
autoplot(fit_fh, type = "estimates")

# Compare two models
fit_sfh <- eblup_sfh(y ~ x1 + x2 + x3, data = mys, vardir = "vardir", W = mys_proxmat)
autoplot(list("Fay-Herriot" = fit_fh, "Spatial FH" = fit_sfh), type = "comparison")

Extract coefficients from a fastsae object

Description

Extract coefficients from a fastsae object

Usage

## S3 method for class 'fastsae'
coef(object, ...)

Arguments

Value

Named vector of estimated coefficients.

Corn and Soybean Survey and Satellite Data in 12 Iowa Counties

Description

Survey and satellite data for corn and soy beans in 12 Iowa counties, originally obtained from the 1978 June Enumerative Survey of the U.S. Department of Agriculture and from LANDSAT satellite observations during the 1978 growing season.

Usage

data(cornsoybean)

Format

A data frame with 37 observations on the following 5 variables:

County

numeric county code.

CornHec

reported hectares of corn from the survey.

SoyBeansHec

reported hectares of soy beans from the survey.

CornPix

number of pixels of corn in the sample segment within county, from satellite data.

SoyBeansPix

number of pixels of soy beans in the sample segment within county, from satellite data.

Details

This dataset is included for demonstration purposes and is originally provided in the sae package.

Source

Battese, G.E., Harter, R.M., and Fuller, W.A. (1988). *An Error-Components Model for Prediction of County Crop Areas Using Survey and Satellite Data.* *Journal of the American Statistical Association*, 83, 28–36.

Corn and Soybean Mean Number of Pixels per Segment for 12 Iowa Counties

Description

County means of number of pixels per segment of corn and soy beans, from satellite data, for 12 counties in Iowa. The dataset includes population size, sample size, and means of auxiliary variables used in the dataset cornsoybean.

Usage

data(cornsoybeanmeans)

Format

A data frame with 12 observations on the following 6 variables:

CountyIndex

numeric county code.

CountyName

name of the county.

SampSegments

number of sample segments in the county (sample size).

PopnSegments

number of population segments in the county (population size).

MeanCornPixPerSeg

mean number of corn pixels per segment in the county.

MeanSoyBeansPixPerSeg

mean number of soy beans pixels per segment in the county.

Details

This dataset is provided for demonstration purposes and is originally distributed with the sae package.

Source

Battese, G.E., Harter, R.M., and Fuller, W.A. (1988). *An Error-Components Model for Prediction of County Crop Areas Using Survey and Satellite Data.* *Journal of the American Statistical Association*, 83, 28–36.

Empirical Best Linear Unbiased Prediction (EBLUP) for the Battese-Harter-Fuller Model

Description

This function estimates small area means or totals using the unit-level model proposed by Battese, Harter, and Fuller (1988), which combines survey data (sample units) and auxiliary population information.

Usage

eblup_bhf(
  formula,
  unit_data,
  Xpop,
  domain_var,
  popsize_var,
  method = c("REML", "ML"),
  popnmean_xpop = NULL,
  B = 100,
  compute_mse = FALSE,
  n_threads = 1,
  seed = -1,
  print_result = TRUE
)

Arguments

Value

List containing EBLUP estimates, fit, and optionally MSE.

References

Battese, G. E., Harter, R. M., and Fuller, W. A. (1988). An error-components model for prediction of county crop areas using survey and satellite data. *Journal of the American Statistical Association*, 83(401), 28-36.

Examples

library(dplyr)
df_meanpop <- cornsoybeanmeans |>
  rename(CornPix = MeanCornPixPerSeg, SoyBeansPix = MeanSoyBeansPixPerSeg)
df_cornsoybean <- cornsoybean |>
  rename(CountyIndex = County)

res <- eblup_bhf(
  formula = CornHec ~ CornPix + SoyBeansPix,
  Xpop = df_meanpop,
  unit_data = df_cornsoybean,
  domain_var = "CountyIndex",
  popsize_var = "PopnSegments"
)

Empirical Best Linear Unbiased Prediction based on a Fay-Herriot Model.

Description

This function gives the Empirical Best Linear Unbiased Prediction (EBLUP) or Empirical Best (EB) predictor under normality based on a Fay-Herriot model.

Usage

eblup_fh(
  formula,
  vardir,
  domain = NULL,
  data,
  method = c("REML", "ML"),
  maxiter = 100,
  precision = 1e-04,
  print_result = TRUE
)

Arguments

Details

The model has a form that is response ~ auxiliary variables. where numeric type response variables can contain NA. When the response variable contains NA it will be estimated with synthetic estimator.

Value

The function returns a list with the following objects: estcoef a data frame with the estimated model coefficients, random_effect_var estimated random effect variance, goodness vector containing several goodness-of-fit measures, df_eblup a data frame that contains y, eblup, random_effect, vardir, mse, and rse.

• y variable response
  

• eblup estimated results for each area
  

• random_effect random effect for each area
  

• vardir variance sampling from the direct estimator for each area
  

• mse Mean Square Error
  

• rse Relative Standart Error (%)
  

References

1. Rao, J. N., & Molina, I. (2015). Small area estimation. John Wiley & Sons.

Examples

library(fastsae)

# Standard Fay-Herriot model
m1 <- eblup_fh(
  y ~ x1 + x2 + x3,
  data = mys,
  vardir = "vardir"
)

Empirical Best Linear Unbiased Prediction based on a Spatial Fay-Herriot Model.

Description

This function gives the Spatial Empirical Best Linear Unbiased Prediction (EBLUP) or Empirical Best (EB) predictor under normality based on a Fay-Herriot model.

Usage

eblup_sfh(
  formula,
  vardir,
  domain = NULL,
  data,
  method = c("REML", "ML"),
  mse_method = c("analytical", "pbmse", "npbmse"),
  W = NULL,
  B = 100,
  n_threads = 1,
  seed = -1,
  maxiter = 100,
  precision = 1e-04,
  print_result = TRUE
)

Arguments

Details

The model has a form that is response ~ auxiliary variables. where numeric type response variables can contain NA. When the response variable contains NA, that domain is treated as unsampled: it is predicted with a full-spatial synthetic (kriging) estimator that borrows strength from sampled neighbors via W, with an analytical MSE approximation.

Value

The function returns a list with the following objects: estcoef a data frame with the estimated model coefficients in the first column (beta), their asymptotic standard errors in the second column (std.error), the t-statistics in the third column (tvalue) and the p-values of the significance of each coefficient in last column (pvalue)
random_effect_var estimated random effect variance
rho estimated spatial autocorrelation parameter
estvarcomp a data frame with parameter, estimate, and std.error
goodness vector containing several goodness-of-fit measures: loglikelihood, AIC, and BIC
df_eblup a data frame that contains the following columns:

• y variable response
  

• eblup estimated results for each area
  

• random_effect random effect for each area
  

• vardir variance sampling from the direct estimator for each area
  

• mse Mean Square Error
  

• rse Relative Standart Error (%)
  

• mse_pb or mse_npb Parametric / Non Parametric bootstrap MSE
  

• mse_pbbc or mse_npbbc Bias Corrected Parametric / Non Parametric bootstrap MSE
  

References

1. Rao, J. N., & Molina, I. (2015). Small area estimation. John Wiley & Sons.

Examples

library(fastsae)

# Spatial Fay-Herriot model
m1 <- eblup_sfh(
  y ~ x1 + x2 + x3,
  data = mys,
  vardir = ~vardir,
  W = mys_proxmat
)

# Spatial Fay-Herriot model with Parametric Bootstrap MSE
m2 <- eblup_sfh(
  y ~ x1 + x2 + x3,
  data = mys,
  vardir = ~vardir,
  mse_method = "pbmse",
  B = 50,
  W = mys_proxmat
)

Empirical Best Linear Unbiased Prediction based on a Spatio-Temporal Fay-Herriot Model.

Description

This function gives the Spatio-Temporal Empirical Best Linear Unbiased Prediction (EBLUP) under normality based on a spatio-temporal Fay-Herriot model. It reimplements the same Fisher-scoring algorithm as eblupSTFH() (package sae, Marhuenda, Molina & Morales 2013), but the estimation loop runs in compiled C++/Armadillo, making it substantially faster and more memory-efficient than the original R implementation – especially for a large number of domains/time periods.

Usage

eblup_stfh(
  formula,
  vardir,
  data,
  domain,
  time,
  W,
  model = c("ST", "S"),
  maxiter = 100,
  precision = 1e-04,
  sigma21_start = NULL,
  rho1_start = 0.5,
  sigma22_start = NULL,
  rho2_start = 0.5,
  compute_mse = FALSE,
  B = 100,
  n_threads = 1,
  seed = -1,
  print_result = TRUE
)

Arguments

Details

This function requires a complete panel (no NA values in the response variable). If your data contains NA values (unsampled areas/time periods), please filter them out before calling this function. Future versions may support automatic handling of unsampled areas.

Value

A list with the same structure as seblup_area(), with additional estvarcomp for spatio-temporal variance/autocorrelation components:

estcoef

data frame with beta, std.error, tvalue, pvalue.

estvarcomp

data frame with estimate, std.error for sigma21, rho1, sigma22, rho2.

goodness

vector with loglike, AIC, BIC.

df_eblup

data frame with eblup, mse, rse, random_effect_u1, random_effect_u2.

When compute_mse = FALSE, mse and rse are NA.

model

model type ("ST" or "S").

convergence

logical, whether the algorithm converged.

n_iter

number of iterations.

B

number of bootstrap replicates (only when compute_mse = TRUE).

References

1. Marhuenda, Y., Molina, I., & Morales, D. (2013). Small area estimation with spatio-temporal Fay-Herriot models. Computational Statistics & Data Analysis, 58, 308-325.

2. Rao, J. N., & Molina, I. (2015). Small area estimation. John Wiley & Sons.

Examples

library(fastsae)
library(dplyr)

mys_panel_nona <- mys_panel |>
  filter(!is.na(y) & year >= 2024)

# Basic EBLUP without MSE
m1 <- eblup_stfh(
  y ~ x1 + x2 + x3,
  data = mys_panel_nona,
  vardir = ~vardir,
  domain = ~area,
  time = ~year,
  W = mys_proxmat[-c(21, 25), -c(21, 25)],
  model = "ST"
)

# EBLUP with parametric bootstrap MSE
m2 <- eblup_stfh(
  y ~ x1 + x2 + x3,
  data = mys_panel_nona,
  vardir = ~vardir,
  domain = ~area,
  time = ~year,
  W = mys_proxmat[-c(21, 25), -c(21, 25)],
  model = "ST",
  compute_mse = TRUE,
  B = 100,
  seed = 42
)

Extract fitted values (EBLUP) from a fastsae object

Description

Extract fitted values (EBLUP) from a fastsae object

Usage

## S3 method for class 'fastsae'
fitted(object, ...)

Arguments

Value

Vector of fitted EBLUP estimates.

mys: mean years of schooling people with disabilities.

Description

A synthetic dataset containing the mean years of schooling for people with disabilities across 42 regencies/municipalities.

Usage

mys

Format

A data frame with 42 rows and 9 variables with 10 domains as non-sampled areas.

area

regency municipality identifier

y

mean years of schooling people with disabilities (NA for unsampled areas)

vardir

variance sampling from the direct estimator for each area

rse

relative standard error (%)

x1

Number of Elementary Schools

x2

Number of Junior High Schools

x3

Number of Senior High Schools

n

Number of eligible samples

weight

Weight

Source

Simulated based on empirical characteristics of the National Socio-Economic Survey (Susenas), BPS-Statistics Indonesia.

mys_panel: mean years of schooling people with disabilities 2016 - 2026.

Description

A synthetic panel dataset containing the mean years of schooling for people with disabilities across 42 regencies/municipalities over 11 time periods (2016 - 2026).

Usage

mys_panel

Format

A data frame with 462 rows and 10 variables.

area

regency municipality identifier

year

year of observation (2016 to 2026)

y

mean years of schooling people with disabilities (NA for unsampled areas)

vardir

variance sampling from the direct estimator for each area

rse

relative standard error (%)

x1

Number of Elementary Schools

x2

Number of Junior High Schools

x3

Number of Senior High Schools

n

Number of eligible samples

weight

Weight

Source

Simulated panel data based on empirical characteristics of the National Socio-Economic Survey (Susenas), BPS-Statistics Indonesia.

Example proximity matrix

Description

A 42 by 42 row-standardized spatial proximity matrix for the 42 areas in mys.

Usage

mys_proxmat

Format

A square numeric matrix with 42 rows and 42 columns with values between 0 and 1.

Source

Simulated based on contiguous administrative boundaries.

Print a fastsae object

Description

Print a fastsae object

Usage

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

Arguments

Value

The original object invisibly.

Print summary of a fastsae object

Description

Print summary of a fastsae object

Usage

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

Arguments

Value

The original object invisibly.

Re-export ggplot2's autoplot generic

Description

These objects are imported from other packages. Follow the links below to see their documentation.

ggplot2

autoplot()

Extract residuals from a fastsae object

Description

Extract residuals from a fastsae object

Usage

## S3 method for class 'fastsae'
residuals(object, ...)

Arguments

Value

Vector of residuals (y - eblup) for sampled areas.

Summarize a fastsae object

Description

Summarize a fastsae object

Usage

## S3 method for class 'fastsae'
summary(object, ...)

Arguments

Value

An object of class summary.fastsae.