R/pr_fitter_nosplit.R

In spencerwoody/saFAB: Selection-adjusted FAB inference

##' Use predictive recursion to obtain a semiparametric estimate of
##' marginal of y (multiple observations per subject, without splitting),
##'
##' @param y sample
##' @param mu0 foobar
##' @param sig0 foobar
##' @param nulltype foobar
##' @return a list, where my_fit is the fitted marginal or y
##'
##' @export
pr_fitter_nosplit <- function(y, mu0=NULL, sig0=NULL, nulltype = "theoretical") {
  ## require(FDRreg)
  ## require(splines)

  ## Split y's
  N <- nrow(y)

  if (!is.matrix(y)) y <- matrix(y, ncol = 1)
  
  yTr <- rowMeans(y) #y's for estimating prior

  ## y's for constructing intervals
  ## yInt <- matrix(y[, -1], nrow = N) #y's for constructing intervals

  yInt <- y

  n <- ncol(yInt)

  ## sampling model of y
  if(nulltype == "theoretical") {
    if(missing(mu0)) mu0 = 0
    if(missing(sig0)) sig0 = 1
  } else if(nulltype == "empirical") {
    null_fit <- efron(yTr, nmids = 200, df = 15,
                     nulltype = "empirical")
    mu0 <- null_fit$mu0
    if(missing(sig0)) sig0 <- null_fit$sig0
  } else {
    error("nulltype must be either \"theoretical\" or \"empirical\"")
  }

  ## fit predictive recursion
  pr_fit <- prfdr(yTr, mu0, sig0)

  ## interpolate on the log scale to get the prior
  logprior_fit <- splinefun(pr_fit$x_grid,
                            log(pr_fit$pitheta_grid),
                            method='natural')
  prior_fit <- function(theta) exp(logprior_fit(theta))

  ## interpolate on the log scale to get the marginal density
  logmy_fit <- splinefun(pr_fit$x_grid, log(pr_fit$fmix_grid), method="natural")

  ## Marginal for data, with no selection
  
  my_fit <- function(y) exp(logmy_fit(y))

  ## Marginal under joint selection

  p <- 1 - pr_fit$pi0

  my_fit_jnt <- function(y, sigma, t, sigma_orig, n = n) {

    int_part <- integrate(int_fun_jnt, 
                          lower = -Inf, 
                          upper = +Inf, n = n,
                          y = y, sigma_orig = sigma_orig,
                          prior_fit = prior_fit, t = t)$value

    (p * int_part + (1 - p) * dnorm(y, 0, sigma_orig / sqrt(n))) *
    ifelse(abs(y) > t, 1, 0)
    
  }

  my_fit_jnt1 <- function(y, sigma, t, sigma_orig, n = n) {

    int_part <- integrate(int_fun_jnt, 
                          lower = -Inf, 
                          upper = +Inf, n = n,
                          y = y, sigma_orig = sigma_orig,
                          prior_fit = prior_fit, t = t)$value

    (p * int_part + (1 - p) * dnorm(y, 0, sigma_orig / sqrt(n))) *
      ifelse(y > t, 1, 0)
    
  }

  ## Marginal under conditional selection
  ## sigma argument needed to match up with Hprime_w_safab function...
  my_fit_cnd <- function(y, sigma, t, sigma_orig, n = n) {
    int_part <- integrate(int_fun_cnd, lower = -Inf, upper = +Inf,
                          y = y, sigma_orig = sigma_orig,
                          prior_fit = prior_fit, n = n, t = t)$value
    
    (p * int_part + 
     (1 - p) * dnorm(y, 0, sigma_orig / sqrt(n)) /
     Pr_S_cnd(0, sigma_orig / sqrt(n), t)) * 
      ifelse(abs(y) > t, 1, 0)
    
  }

  my_fit_cnd1 <- function(y, sigma, t, sigma_orig, n = n) {
    int_part <- integrate(int_fun_cnd1, lower = -Inf, upper = +Inf,
                          y = y, sigma_orig = sigma_orig,
                          prior_fit = prior_fit, n = n, t = t)$value
    
    (p * int_part + 
     (1 - p) * dnorm(y, 0, sigma_orig / sqrt(n)) /
     Pr_S_cnd1(0, sigma_orig / sqrt(n), t)) * 
      ifelse(y > t, 1, 0)
    
  }

  return(list(
    prior_fit = prior_fit,
    my_fit = my_fit,
    my_fit_jnt = my_fit_jnt,
    my_fit_jnt1 = my_fit_jnt1,
    my_fit_cnd = my_fit_cnd,
    my_fit_cnd1 = my_fit_cnd1,
    pr_fit = pr_fit,
    mu0 = mu0,
    sig0 = sig0,
    phat = p
  ))

}