Title:	Gaussian Processes for Social Science
Version:	1.0.3
Date:	2026-04-01
Description:	Provides Gaussian process (GP) regression tools for social science inference problems. GPs combine flexible nonparametric regression with principled uncertainty quantification: rather than committing to a single model fit, the posterior reflects lesser knowledge at the edge of or beyond the observed data, where other approaches become highly model-dependent. The package reduces user-chosen hyperparameters from three to zero and supplies convenience functions for regression discontinuity (gp_rdd()), interrupted time-series (gp_its()), and general GP fitting (gpss(), gp_train(), gp_predict()). Methods are described in Cho, Kim, and Hazlett (2026) <doi:10.1017/pan.2026.10032>.
License:	GPL (≥ 3)
Encoding:	UTF-8
RoxygenNote:	7.3.3
URL:	https://doeunkim.org/gpss/, https://github.com/doeun-kim/gpss
BugReports:	https://github.com/doeun-kim/gpss/issues
Imports:	Rcpp (≥ 1.0.14), Matrix, ggplot2, rlang
LinkingTo:	Rcpp, RcppArmadillo
Depends:	R (≥ 3.5.0)
LazyData:	true
Suggests:	MASS, posterior, testthat (≥ 3.0.0)
NeedsCompilation:	yes
Packaged:	2026-04-01 20:02:06 UTC; chadhazlett
Author:	Soonhong Cho [aut], Doeun Kim [aut], Chad Hazlett [aut, cre]
Maintainer:	Chad Hazlett <chazlett@ucla.edu>
Repository:	CRAN
Date/Publication:	2026-04-08 19:10:17 UTC

gpss: Gaussian Processes for Social Science

Description

Provides Gaussian process (GP) regression tools for social science inference problems. GPs combine flexible nonparametric regression with principled uncertainty quantification: rather than committing to a single model fit, the posterior reflects lesser knowledge at the edge of or beyond the observed data, where other approaches become highly model-dependent. The package reduces user-chosen hyperparameters from three to zero and supplies convenience functions for regression discontinuity (gp_rdd()), interrupted time-series (gp_its()), and general GP fitting (gpss(), gp_train(), gp_predict()). Methods are described in Cho, Kim, and Hazlett (2026) doi:10.1017/pan.2026.10032.

Usage

gpss(
  formula,
  data,
  b = NULL,
  s2 = 0.3,
  optimize = FALSE,
  scale = TRUE,
  prior_mean = NULL
)

Arguments

Value

Author(s)

Maintainer: Chad Hazlett chazlett@ucla.edu

Authors:

• Soonhong Cho tnsehdtm@gmail.com

• Doeun Kim doeun2@ucla.edu (ORCID)

See Also

Useful links:

• https://doeunkim.org/gpss/

• https://github.com/doeun-kim/gpss

• Report bugs at https://github.com/doeun-kim/gpss/issues

Examples

# -- Basic fitting and prediction -----------------------------
library(gpss)
data(lalonde)
dat <- transform(lalonde, race_ethnicity = factor(race_ethnicity))

idx <- sample(seq_len(nrow(dat)), 500)
mod <- gpss(re78 ~ nsw + age + educ + race_ethnicity, data = dat[idx, ])
p   <- predict(mod, dat[-idx, ])
head(p)

# -- G-computation for the ATT (LaLonde data) ----------------
# The LaLonde dataset pairs 185 treated units from the NSW
# job-training program with 2490 PSID control units.  The
# experimental benchmark ATT is about $1794, but the naive
# difference in means is severely biased (~-$15k) due to
# poor covariate overlap.
#
# Strategy: fit a GP on the control group to learn E[Y(0)|X],
# predict counterfactual Y(0) for the treated, and compute
# ATT = mean(Y_obs - Y0_hat) with a Neyman-style SE.
# See Cho, Kim, and Hazlett (2026, Political Analysis).

# Use gp_train / gp_predict for direct control over covariates
cat_vars <- c("black", "hisp", "married", "nodegr", "u74", "u75")
all_vars <- c("age", "educ", "re74", "re75",
              "black", "hisp", "married", "nodegr", "u74", "u75")

X <- lalonde[, all_vars]
Y <- lalonde$re78
D <- lalonde$nsw

# Fit GP on control group (optimize s2 via MLE)
gp_mod <- gp_train(X = X[D == 0, ], Y = Y[D == 0],
                   optimize = TRUE,
                   mixed_data = TRUE, cat_columns = cat_vars)

# Predict E[Y(0) | X] for treated units
gp_pred <- gp_predict(gp_mod, Xtest = X[D == 1, ])

Y1     <- Y[D == 1]                    # observed Y(1)
Y0_hat <- gp_pred$Ys_mean_orig         # predicted E[Y(0)|X]
tau_i  <- Y1 - Y0_hat                  # unit-level effects
att    <- mean(tau_i)                   # ATT point estimate

# -- Neyman-style SE for the ATT --
# The estimator is a difference of two means:
# ATT = mean(Y_obs - Y0_hat).
# So a Neyman-style SE:
# Var(mean(Y_obs)) + Var(mean(Y0_hat)).
#
# Two sources of uncertainty:
#
# 1. Sampling variance of observed outcomes: var(Y1) / n1
#
# 2. Posterior uncertainty in the counterfactuals: because Y(0)
#    predictions are correlated across units through the GP
#    posterior, we need the full covariance matrix V0 = f_cov_orig
#    (the posterior covariance of E[Y(0)|X] at treated X values).
#    Its contribution to Var(ATT) is (1/n1^2) * sum(V0),
#    i.e. 1'V0 1 / n1^2.

n1     <- sum(D)
V0     <- gp_pred$f_cov_orig
se_att <- sqrt(var(Y1) / n1 + sum(V0) / n1^2)

cat(sprintf("GP ATT:  $%.0f (SE = $%.0f)\n", att, se_att))
cat(sprintf("95%% CI: [$%.0f, $%.0f]\n",
            att - 1.96 * se_att, att + 1.96 * se_att))
# GP estimate is close to the experimental benchmark (~$1794);
# the wide CI reflects genuine extrapolation uncertainty.

Internal function for placebo checks

Description

Internal function for placebo checks

Usage

.run_placebo_checks(
  y,
  X_full,
  n_pre,
  kernel_type,
  period,
  time_col,
  scale,
  interval_type,
  alpha,
  optimize,
  mixed_data,
  cat_columns,
  placebo_periods,
  verbose
)

Gaussian Process Interrupted Time Series

Description

Gaussian Process Interrupted Time Series

Usage

gp_its(
  y,
  dates,
  date_treat,
  covariates = NULL,
  kernel_type = "gaussian",
  b = NULL,
  s2 = 0.3,
  period = NULL,
  auto_period = FALSE,
  time_col = NULL,
  scale = TRUE,
  optimize = TRUE,
  interval_type = c("prediction", "confidence"),
  alpha = 0.05,
  mixed_data = FALSE,
  cat_columns = NULL,
  placebo_check = FALSE,
  placebo_periods = NULL,
  verbose = FALSE
)

Arguments

Value

An object of class "gp_its", a list containing:

Examples

# Simulated interrupted time series
set.seed(42)
n <- 50
dates <- seq(as.Date("2020-01-01"), by = "month", length.out = n)
y <- rnorm(n) + c(rep(0, 30), rep(2, 20))
res <- gp_its(y, dates, date_treat = as.Date("2022-07-01"))
print(res)

gp_predict

Description

to predict outcome values at testing points by feeding the results obtained from gp_train()

Usage

gp_predict(gp, Xtest, prior_mean = NULL)

Arguments

Value

Examples

data(lalonde)
cat_vars <- c("race_ethnicity", "married")
all_vars <- c("age","educ","re74","re75","married", "race_ethnicity")

X <- lalonde[,all_vars]
Y <- lalonde$re78
D <- lalonde$nsw

X_train <- X[D==0,]
Y_train <- Y[D==0]
X_test <- X[D == 1,]
Y_test <- Y[D == 1]

gp_train.out <- gp_train(X = X_train, Y = Y_train,
optimize=TRUE, mixed_data = TRUE,
cat_columns = cat_vars)
gp_predict.out <- gp_predict(gp_train.out, X_test)

gp_rdd

Description

This function implements a regression discontinuity (RD) design using Gaussian Process (GP) regression on each side of the cutoff. Separate GP models are estimated for observations below and above the threshold, allowing for flexible and fully nonparametric functional forms. The treatment effect is computed as the difference between the predicted conditional means at the cutoff from the right and left limits. Standard errors and confidence intervals are constructed using the posterior predictive variance from the two independently fitted GPs.

Usage

gp_rdd(
  X,
  Y,
  cut,
  alpha = 0.05,
  b = NULL,
  trim = FALSE,
  trim_k_value = 0.1,
  scale = TRUE
)

Arguments

Value

Examples

n <- 100
tau <- 3
cut <- 0
x <- rnorm(n, 0, 1)
y <- rnorm(n, 0, 1) + tau*ifelse(x>cut, 1, 0)
gp_rdd.out <- gp_rdd(x, y, cut)
gp_rdd_plot(gp_rdd.out)

gp_rdd_plot

Description

to draw an RD plot using the results obtained from gp_rdd()

Usage

gp_rdd_plot(gp_rdd_res, l_col = "blue", r_col = "red")

Arguments

Value

A ggplot object showing the RD plot.

Examples

library(ggplot2)
n <- 100
tau <- 3
cut <- 0
x <- rnorm(n, 0, 1)
y <- rnorm(n, 0, 1) + tau*ifelse(x>cut, 1, 0)
gp_rdd.out <- gp_rdd(x, y, cut)
gp_rdd_plot(gp_rdd.out) +
 geom_vline(xintercept = cut, lty="dashed")

gp_train

Description

to train GP model with training data set

Usage

gp_train(
  X,
  Y,
  b = NULL,
  s2 = 0.3,
  optimize = FALSE,
  scale = TRUE,
  kernel_type = "gaussian",
  period = NULL,
  time_col = NULL,
  mixed_data = FALSE,
  cat_columns = NULL,
  Xtest = NULL,
  prior_mean = NULL
)

Arguments

Value

Examples

data(lalonde)
cat_vars <- c("race_ethnicity", "married")
all_vars <- c("age","educ","re74","re75","married", "race_ethnicity")

X <- lalonde[,all_vars]
Y <- lalonde$re78
D <- lalonde$nsw

X_train <- X[D==0,]
Y_train <- Y[D==0]
X_test <- X[D == 1,]
Y_test <- Y[D == 1]

gp_train.out <- gp_train(X = X_train, Y = Y_train,
optimize=TRUE, mixed_data = TRUE,
cat_columns = cat_vars)
gp_predict.out <- gp_predict(gp_train.out, X_test)

Data from National Supported Work program and Panel Study in Income Dynamics

Description

Dehejia and Wahba (1999) sample of data from Lalonde (1986). This data set includes 185 treated units from the National Supported Work (NSW) program, paired with 2490 control units drawn from the Panel Study of Income Dynamics (PSID-1).

The treatment variable of interest is nsw, which indicates that an individual was in the job training program. The main outcome of interest is real earnings in 1978 (re78). The remaining variables are characteristics of the individuals, to be used as controls.

Usage

lalonde

Format

A data frame with 2675 rows and 14 columns.

nsw

treatment indicator: participation in the National Supported Work program.

re78

real earnings in 1978 (outcome)

u78

unemployed in 1978; actually an indicator for zero income in 1978

age

age in years

black

indicator for identifying as black

hisp

indicator for identifying as Hispanic

race_ethnicity

factor for self-identified race/ethnicity; same information as black and hisp in character form.

married

indicator for being married

re74

real income in 1974

re75

real income in 1975

u74

unemployment in 1974; actually an indicator for zero income in 1974

u75

unemployment in 1975; actually an indicator for zero income in 1975

educ

Years of education of the individual

nodegr

indicator for no high school degree; actually an indicator for years of education less than 12

References

Dehejia, Rajeev H., and Sadek Wahba. "Causal effects in non-experimental studies: Reevaluating the evaluation of training programs." Journal of the American statistical Association 94.448 (1999): 1053-1062.

LaLonde, Robert J. "Evaluating the econometric evaluations of training programs with experimental data." The American economic review (1986): 604-620.

predict method for gpss objects

Description

predict method for gpss objects

Usage

## S3 method for class 'gpss'
predict(
  object,
  newdata = NULL,
  type = "response",
  format = "default",
  interval = "confidence",
  level = 0.95,
  prior_mean = NULL,
  ...
)

Arguments

Value

If format = "default", a numeric matrix with columns fit, lwr, and upr giving the posterior mean and the lower and upper bounds of the requested interval for each row of newdata. If format = "rvar", a copy of newdata with an additional column rvar containing a posterior random variable (class rvar from the posterior package).

Examples

library(gpss)
data(lalonde)
# categorical variables must be encoded as factors
dat <- transform(lalonde, race_ethnicity = factor(race_ethnicity))
# train and test sets
idx <- sample(seq_len(nrow(dat)), 500)
dat_train <- dat[idx, ]
dat_test <- dat[-idx, ]
# Fit model
mod <- gpss(re78 ~ nsw + age + educ + race_ethnicity, data = dat_train)
p <- predict(mod, newdata = dat_test)
p_confidence99 <- predict(mod, newdata = dat_test, interval = "confidence", level = 0.99)

Print method for gp_its objects

Description

Print method for gp_its objects

Usage

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

Arguments

Value

The input object x, returned invisibly.

Print method for gp_rdd objects

Description

Print method for gp_rdd objects

Usage

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

Arguments

Value

The input object x, returned invisibly.

Summary method for gp_its objects

Description

Summary method for gp_its objects

Usage

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

Arguments

Value

The input object object, returned invisibly.

Summary method for gp_rdd objects

Description

Summary method for gp_rdd objects

Usage

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

Arguments

Value

The input object object, returned invisibly.