# R/p-value_calculation.R In davidruegamer/coinflibs: Conditional Inference after Likelihood-based Selection

```mult_tnorm_surv <- function(x, mean, sd, limits, twosided = TRUE)
{

limits <- limits[order(limits[,1]),]
pnormstd <- function(x) pnorm((x - mean)/sd)
limitX <- which(limits[,2] > x & limits[,1] < x)

# allow for multiple means
if(length(mean)==1)
{

denom <-
sum(sapply(1:nrow(limits), function(i)
pnormstd(limits[i,2]) - pnormstd(limits[i,1])))

}else{

denom <- sapply(1:nrow(limits), function(i)
pnormstd(limits[i,2]) - pnormstd(limits[i,1]))
denom <- apply(denom, 1, sum)

}

limitXb <- limits[limitX,2]

if(length(mean)==1)
{

massLargerXb <- if(nrow(limits)>limitX)
sum(sapply((limitX+1):nrow(limits), function(i)
pnormstd(limits[i,2]) - pnormstd(limits[i,1]))) else 0

}else{

if(nrow(limits)>limitX)
{

massLargerXb <- sapply((limitX+1):nrow(limits), function(i)
pnormstd(limits[i,2]) - pnormstd(limits[i,1]))
massLargerXb <- apply(massLargerXb, 1, sum)

}else{

massLargerXb <- rep(0, length(mean))

}

}

nom <- (massLargerXb + pnormstd(limitXb)) - pnormstd(x)

this_pval <- nom / denom
this_pval[is.nan(this_pval) & mean < x] <- 0
this_pval[is.nan(this_pval) & mean > x] <- 1

if(all(is.na(this_pval) | is.nan(this_pval)))
this_pval <- bryc.tnorm.surv(z = x, sd = sd,
a = limits[limitX,1],
b = limits[limitX,2])

if(twosided) this_pval <- sapply(this_pval,
function(x) 2*min(c(x, 1-x)))

return(this_pval)

}

#' Calculate confidence interval for truncated normal variable
#'
#' @param meanEst mean
#' @param sd standard deviation
#' @param limits \code{Intervals} object containing all limits
#' @param alpha value specifying the interval coverage, which is \code{(1-alpha)}
#' @param gridpts,griddepth parameters for grid search used in the confidence interval
#' function of the package \code{selectiveInference}.
#' @param range gives the range \code{range*sd}, in which the quantile values are searched
#'
#' @examples
#' library(intervals)
#' # one restriction -- in this case equal to a call to \code{\link{msm::qtnorm}}
#' int1 <- Intervals(c(-3,1))
#' mean = -0.2
#' sd = 1
#' ci_tnorm(meanEst = mean, sd = sd, limits = int1, alpha = 0.05)
#'
#' # two intervals -- 1-alpha coverage not possible
#' (int2 <- Intervals(matrix(c(-3,-2,-1,1),
#'                           ncol = 2, byrow = TRUE)))
#' ci_tnorm(meanEst = mean, sd = sd, limits = int2, alpha = 0.05)
#'
#' # two intervals -- 1-alpha coverage possible, but interval not symmetric
#' (int3 <- Intervals(matrix(c(-3,-2.5,-1.9,1),
#'                           ncol = 2, byrow =  TRUE)))
#' ci_tnorm(meanEst = mean, sd = sd, limits = int3, alpha = 0.05)
#'
#' @importFrom msm qtnorm
#' @export
#'
ci_tnorm <- function(meanEst, sd, limits, alpha,
gridpts = 500, griddepth = 3, range = 100)
{

# grid search adapted from the selectiveInference package

limitsX <- which(limits[,2] > meanEst & limits[,1] < meanEst)

xg = seq(-range * sd, range * sd, length = gridpts)

if(nrow(limits)==1){

fun = function(x) {
tnorm.surv(z = meanEst, mean = x, sd = sd,
a = limits[limitsX, 1],
b = limits[limitsX, 2]) }

}else{

fun = function(x) {
mult_tnorm_surv(x = meanEst, mean = x, sd = sd,
limits = limits, twosided = FALSE) }

}

int = grid.search(xg, fun, alpha/2, 1-alpha/2,
gridpts, griddepth)

return(Intervals(int))

}
```
davidruegamer/coinflibs documentation built on May 14, 2019, 10:33 a.m.