ebTobit {ebTobit}

Empirical Bayes Matrix Estimation under a Tobit Likelihood

Description

Fit and estimate the nonparametric maximum likelihood estimator in R^p (p >= 1) when the likelihood is Gaussian and possibly interval censored. If p = 1, then L, R, and gr may be vectors (they are immediately converted into matrices internally).

Usage

ebTobit(
  L,
  R = L,
  gr = (R + L)/2,
  s1 = 1,
  algorithm = "EM",
  pos_lik = TRUE,
  ...
)

Arguments

L

n x p matrix of lower bounds on observations

R

n x p matrix of upper bounds on observations

gr

m x p matrix of grid points

s1

a single numeric standard deviation or an n x p matrix of standard deviations

algorithm

method to fit prior, either a function or function name

pos_lik

boolean indicating whether to lower-bound the likelihood matrix with .Machine$double.xmin (default: TRUE); helps avoid possible divide-by-zero errors in algorithm

...

further arguments passed into fitting method such as rtol and maxiter, see for example EM

Details

Each observation is stored in a pair of matrices, L and R. If L_ij = R_ij then a direct measurement X_ij ~ N(theta, s1^2) is made; if L_ij < R_ij then the measurement is censored so that L_ij < X_ij < R_ij.

To use a custom fitting algorithm, define a function MyAlg that takes in an n x m likelihood matrix: P_ij = P(L_i, R_i | theta = t_j) and returns a vector of estimated prior weights for t_j. Once MyAlg is defined, fit the prior by using algorithm = "MyAlg" or use the function itself algorithm = MyAlg.

Alternative fitting algorithms "ConvexPrimal"and "ConvexDual" have been (wrappers of REBayes::KWPrimal and REBayes::KWDual, respectively) included and can be used if MOSEK and REBayes are properly installed.

Value

a fitted ebTobit object containing at least the prior weights, corresponding grid/support points, and likelihood matrix relating the grid to the observations

Examples

set.seed(1)
n <- 100
p <- 5
r <- 2
U.true <- matrix(stats::rexp(n*r), n, r)
V.true <- matrix(sample(x = c(1,4,7), 
                        size = p*r, 
                        replace = TRUE, 
                        prob = c(0.7, 0.2, 0.1)), 
                 p, r)
TH <- tcrossprod(U.true, V.true)
X <- TH + matrix(stats::rnorm(n*p), n, p)

# fit uncensored method
fit1 <- ebTobit(X)

# fit left-censored method
ldl <- 1 # lower and upper detection limits
udl <- Inf
L <- ifelse(X < ldl, 0, ifelse(X <= udl, X, udl))
R <- ifelse(X < ldl, ldl, ifelse(X <= udl, X, Inf))
fit2 <- ebTobit(L, R)