Package 'fastsae'

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

object	An object of class fastsae, or a (named) list of fastsae objects.
type	Type of plot to create. For a single fastsae object: "comparison", "mse", "estimates", or "scatter". For a list of fastsae objects: "comparison", "mse", or "scatter".
...	Additional arguments passed to internal plotting helpers (e.g. title) or to ggplot2 layers.

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

object	An object of class fastsae.
...	Additional 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

formula	An object of class 'formula' describing the model.
unit_data	A 'data.frame' containing the unit-level survey data.
Xpop	A 'data.frame' containing auxiliary variables and domain info.
domain_var	A character string giving the column name for domain identifier.
popsize_var	A character string for population size variable.
method	Fitting method: "REML" (default) or "ML".
popnmean_xpop	Population mean of auxiliary variables per domain.
B	Number of bootstrap replicates for MSE (if compute_mse = TRUE).
compute_mse	If TRUE, compute bootstrap MSE.
n_threads	Number of threads for parallel computation.
seed	Random seed for reproducibility.
print_result	Print results (default TRUE).

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

formula	an object of class formula that contains a description of the model to be fitted.
vardir	vector or column names from data that contain variance sampling from the direct estimator.
domain	vector, column name or one-sided formula referencing a domain names column in data. If NULL, the domains are numbered consecutively.
data	a data frame or a data frame extension (e.g. a tibble).
method	Fitting method can be chosen between 'ML' and 'REML'.
maxiter	maximum number of iterations allowed in the Fisher-scoring algorithm.
precision	convergence tolerance limit for the Fisher-scoring algorithm.
print_result	print coefficient or not, default value is TRUE.

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

formula	an object of class formula that contains a description of the model to be fitted. The variables included in the formula must be contained in the data.
vardir	vector or column names from data that contain variance sampling from the direct estimator for each area.
domain	vector, column name or one-sided formula referencing a domain names column in data. If NULL, the domains are numbered consecutively.
data	a data frame or a data frame extension (e.g. a tibble).
method	Fitting method can be chosen between 'ML' and 'REML'.
mse_method	a character string determining the estimation method of the MSE. Methods that can be chosen: "analytical", "pbmse", or "npbmse". When there are unsampled domains, "pbmse"/"npbmse" bootstrap MSE is only defined for the sampled domains; unsampled domains automatically get a full-spatial synthetic (kriging) prediction and analytical MSE merged into the result regardless of mse_method.
W	A square matrix with dimension equal to the TOTAL number of domains in data (including any unsampled domains where the response is NA). It should contain the row-standardized spatial weights (proximities) between ALL domains, with values ranging from 0 to 1. Rows and columns must be ordered consistently with the domain identifiers in data. Do NOT pre-subset W to sampled domains only – unsampled domains' spatial MSE/ prediction relies on their relation (in W) to sampled neighbors.
B	Number of bootstrap replications when mse_method = "pbmse" or "npbmse".
n_threads	Number of threads used in parallel computation (default 1). Values less than or equal to 0 use all available cores.
seed	Integer seed for bootstrap resampling. A value of -1 leaves the current R RNG state unchanged.
maxiter	maximum number of iterations allowed in the Fisher-scoring algorithm. Default is 100 iterations.
precision	convergence tolerance limit for the Fisher-scoring algorithm. Default value is 0.0001.
print_result	print coefficient or not, default value is TRUE.

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

formula	an object of class formula describing the model to fit (response ~ auxiliary variables). Variables must be present in data.
vardir	vector, column name or one-sided formula referencing a column in data, with the sampling variances of the direct estimator.
data	a data frame (or extension) with domain * time rows, sorted so that all time periods of domain 1 come first, then all periods of domain 2, and so on (i.e. domain-major order) – exactly as required by eblupSTFH().
domain	vector, column name or one-sided formula referencing a domain names column in data. If NULL, the domains are numbered consecutively.
time	vector, column name, or one-sided formula referencing a time names column in data.
W	a square proximity/spatial weights matrix of dimension domain x domain (row-standardized, values typically in $[0,1]$).
model	character, either "ST" (spatio-temporal, default) or "S" (spatial only, no AR(1) temporal component).
maxiter	maximum number of Fisher-scoring iterations. Default 100.
precision	convergence tolerance for the Fisher-scoring algorithm. Default 1e-4.
sigma21_start, rho1_start, sigma22_start, rho2_start	starting values for the variance/autocorrelation components. Defaults mirror eblupSTFH(): 0.5 * median(vardir) for the variances and 0.5 for the autocorrelations. rho2_start is ignored when model = "S".
compute_mse	logical, if TRUE computes parametric bootstrap MSE using B bootstrap replicates. Default FALSE.
B	number of bootstrap replicates for MSE computation. Only used when compute_mse = TRUE. Default 100.
n_threads	number of threads for parallel bootstrap MSE computation. Use 0 for all available cores. Default 1.
seed	random seed for bootstrap. Use -1 for no seed. Default -1.
print_result	print the estimated coefficients or not. Default TRUE.

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

object	An object of class fastsae.
...	Additional 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

x	An object of class fastsae.
...	Additional arguments passed to print methods.

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

x	An object of class summary.fastsae.
...	Additional arguments.

Value

The original object invisibly.

---

Extract residuals from a fastsae object

Description

Extract residuals from a fastsae object

Usage

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

Arguments

object	An object of class fastsae.
...	Additional 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

object	An object of class fastsae.
...	Additional arguments.

Value

An object of class summary.fastsae.

---