twostageTE (version 1.3)

threshold_estimate_ir: Threshold estimate based on IR

Description

Uses isotonic regression and PAVA to form a point estimate.

Usage

threshold_estimate_ir(explanatory, response, Y_0)

Arguments

explanatory

Explanatory sample points

response

Observed responses at the explanatory sample points

Y_0

Threshold of interest

Value

list(threshold_estimate_explanatory = estim_x, threshold_estimate_response = fit$y[ind], threshold = Y_0, Y_hat = fit$y, index = ind)

threshold_estimate_explanatory

Point estimate of d_0

threshold_estimate_response

Estimate of f(d_0), which may not be exactly equal to the desired threshold

threshold

Threshold of interest (equal to Y_0 input)

Y_hat

Fitted values from PAVA

index

index that corresponds to the point estimate, so that Y_hat[index]=threshold_estimate_response

Details

This is an internal function not meant to be called directly. It function relies on the PAVA algorithm to form a point estimate.

Examples

X=runif(25, 0,1)
Y=X^2+rnorm(n=length(X), sd=0.1)
stageOneAnalysis(X, Y, 0.25, type="IR-wald", 0.99)

## The function is currently defined as
function (explanatory, response, Y_0) 
{
    n = length(response)
    if (sum(response < Y_0) == n) {
        warning("Y_0 is outside observed region")
        list(threshold_estimate_explanatory = max(explanatory), 
            threshold_estimate_response = max(response), threshold = Y_0, 
            Y_hat = max(response), index = n)
    }
    else if (sum(response >= Y_0) == n) {
        warning("Y_0 is outside observed region")
        list(threshold_estimate_explanatory = min(explanatory), 
            threshold_estimate_response = min(response), threshold = Y_0, 
            Y_hat = min(response), index = 1)
    }
    else {
        fit = pava(explanatory, response)
        if (sum(fit$y >= Y_0) == 0) {
            warning("estimate is on the boundary")
            ind = n
            estim_x = fit$x[ind]
        }
        else if (sum(fit$y <= Y_0) == 0) {
            warning("estimate is on the boundary")
            ind = min(which(fit$y >= Y_0))
            estim_x = fit$x[ind]
        }
        else {
            ind = min(which(fit$y >= Y_0))
            estim_x = fit$x[ind]
        }
        list(threshold_estimate_explanatory = estim_x, 
            threshold_estimate_response = fit$y[ind], 
            threshold = Y_0, Y_hat = fit$y, index = ind)
    }
  }