hierSDR (version 0.1)

hier.sphd: Main hierarchical sufficient dimension reduction fitting function

Description

fits hierarchically nested sufficient dimension reduction models

Usage

hier.sphd(
  x,
  y,
  z,
  z.combinations,
  d,
  weights = rep(1L, NROW(y)),
  maxit = 250L,
  tol = 1e-09,
  h = NULL,
  opt.method = c("lbfgs2", "lbfgs.x", "bfgs.x", "bfgs", "lbfgs", "spg", "ucminf", "CG",
    "nlm", "nlminb", "newuoa"),
  init.method = c("random", "phd"),
  vic = TRUE,
  grassmann = TRUE,
  nn = NULL,
  nn.try = c(0.15, 0.25, 0.5, 0.75, 0.9, 0.95),
  n.random = 100L,
  optimize.nn = FALSE,
  separate.nn = FALSE,
  constrain.none.subpop = TRUE,
  verbose = TRUE,
  degree = 2,
  pooled = FALSE,
  maxk = 5000,
  ...
)

Arguments

x

an n x p matrix of covariates, where each row is an observation and each column is a predictor

y

vector of responses of length n

z

an n x C matrix of binary indicators, where each column is a binary variable indicating the presence of a binary variable which acts as a stratifying variable. Each combination of all columns of z pertains to a different subpopulation. WARNING: do not use too many binary variables in z or else it will quickly result in subpopulations with no observations

z.combinations

a matrix of dimensions 2^C x C with each row indicating a different combination of the possible values in z. Each combination represents a subpopulation. This is necessary because we need to specify a different structural dimension for each subpopulation, so we need to know the ordering of the subpopulations so we can assign each one a structural dimension

d

an integer vector of length 2^C of structural dimensions. Specified in the same order as the rows in z.combinations

weights

vector of observation weights

maxit

maximum number of iterations for optimization routines

tol

convergence tolerance for optimization routines. Defaults to 1e-6

h

bandwidth parameter. By default, a reasonable choice is selected automatically

opt.method

optimization method to use. Available choices are c("lbfgs2", "lbfgs.x", "bfgs.x", "bfgs", "lbfgs", "spg", "ucminf", "CG", "nlm", "nlminb", "newuoa")

init.method

method for parameter initialization. Either "random" for random initialization or "phd" for a principle Hessian directions initialization approach

vic

logical value of whether or not to compute the VIC criterion for dimension determination

grassmann

logical value of whether or not to enforce parameters to be on the Grassmann manifold

nn

nearest neighbor parameter for locfit.raw

nn.try

vector of nearest neighbor parameters for locfit.raw to try in random initialization

n.random

integer number of random initializations for parameters to try

optimize.nn

should nn be optimized? Not recommended

separate.nn

should each subpopulation have its own nn? If TRUE, optimization takes much longer. It is rarely better, so recommended to set to FALSE

constrain.none.subpop

should the "none" subpopulation be constrained to be contained in every other subpopulation's dimension reduction subspace? Recommended to set to TRUE

verbose

should results be printed along the way?

degree

degree of kernel to use

pooled

should the estimator be a pooled estimator?

maxk

maxk parameter for locfit.raw. Set to a large number if an out of vertex space error occurs.

...

extra arguments passed to locfit.raw

Value

A list with the following elements

• beta a list of estimated sufficient dimension reduction matrices, one for each subpopulation

• beta.init a list of the initial sufficient dimension reduction matrices, one for each subpopulation -- do not use, just for the sake of comparisons

• directions a list of estimated sufficient dimension reduction directions (i.e. the reduced dimension predictors/variables), one for each subpopulation. These have number of rows equal to the sample size for the subpopulation and number of columns equal to the specified dimensions of the reduced dimension spaces.

• y.list a list of vectors of responses for each subpopulation

• z.combinations the z.combinations specified as an input

• cov list of variance covariance matrices for the covariates for each subpopulation

• sqrt.inv.cov list of inverse square roots of the variance covariance matrices for the covariates for each subpopulation. These are used for scaling

• solver.obj object returned by the solver/optimization function

• value value of the objective function at the solution

• value.init value of the objective function at the initial beta (beta.init) used

• vic.est.eqn the average (unpenalized) VIC value across the r different input values. This assesses model fit

• vic.eqns the individual (unpenalized) VIC values across the r input values. Not used.

• vic the penalized VIC value. This is used for dimension selection, with dimensions chosen by the set of dimensions that minimize this penalized vic value that trades off model complexity and model fit

Examples

# NOT RUN {
library(hierSDR)

set.seed(123)
dat <- simulate_data(nobs = 200, nvars = 6,
                     x.type = "some_categorical",
                     sd.y = 1, model = 2)

x <- dat$x ## covariates
z <- dat$z ## factor indicators
y <- dat$y ## response

dat$beta ## true coefficients that generate the subspaces

dat$z.combinations ## what combinations of z represent different subpops

## correct structural dimensions:
dat$d.correct

## fit hier SPHD model:

# }
# NOT RUN {
hiermod <- hier.sphd(x, y, z, dat$z.combinations, d = dat$d.correct,
                     verbose = FALSE, maxit = 250, maxk = 8200)

## validated inf criterion for choosing dimensions (the smaller the better)
hiermod$vic


cbind(hiermod$beta[[4]], NA, dat$beta[[4]])

## angles between estimated and true subspaces for each population:
mapply(function(x,y) angle(x,y), hiermod$beta, dat$beta)

## projection difference norm between estimated and true subspaces for each population:
mapply(function(x,y) projnorm(x,y), hiermod$beta, dat$beta)
# }
# NOT RUN {

# }