Skip to content

Commit 95e5e00

Browse files
timothyfrasertimothyfraser
authored
Create functions_k.R
1 parent ad21b1b commit 95e5e00

1 file changed

Lines changed: 164 additions & 0 deletions

File tree

‎functions/functions_k.R‎

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
#' @name qk
2+
#' @title k-factor Quantiles
3+
#' @description
4+
#' Function to return quantiles for k-factors.
5+
#' Intended for estimating confidence intervals for failure rates.
6+
#' @param p:[dbl] vector of probabilities / percentile(s)
7+
#' @param r:int number of failures (non-negative integers; can include zero)
8+
#' @param .time:logical logical; is case time-censored data?
9+
#' @param .failure:logical logical; is case failure-censored data?
10+
#'
11+
#' @importFrom dplyr `case_when`
12+
qk = function(p, r, .time = FALSE, .failure = FALSE){
13+
# Testing values
14+
# p = 0.95; r = 20; .time = FALSE; .failure = FALSE
15+
16+
# Input error handling
17+
stopifnot(is.numeric(p) & p >= 0 & p <= 1)
18+
stopifnot(is.numeric(as.integer(r)))
19+
stopifnot(is.logical(.time))
20+
stopifnot(is.logical(.failure))
21+
stopifnot(
22+
(.time == TRUE & .failure == FALSE) |
23+
(.time == FALSE & .failure == FALSE) |
24+
(.time == FALSE & .failure == TRUE) )
25+
stopifnot(
26+
# R should either be greater than 0
27+
(r > 0) |
28+
# OR
29+
# zero and the data should be time censored
30+
(r == 0 & .time == TRUE)
31+
)
32+
33+
# Evaluate if p is in the upper or lower tail
34+
.upper = p > 0.5
35+
36+
# Does r == 0?
37+
.zerofailures = r == 0
38+
39+
40+
k = case_when(
41+
# 1+ failures AND complete data AND UPPER tail --> Get k-factor for r as normal
42+
.zerofailures == FALSE & .time == FALSE & .failure == FALSE & .upper == TRUE ~ qchisq(p, df = 2*r) / (2*r),
43+
# 1+ failures AND complete data AND LOWER tail --> Get k-factor for r as normal
44+
.zerofailures == FALSE & .time == FALSE & .failure == FALSE & .upper == FALSE ~ qchisq(p, df = 2*r) / (2*r),
45+
# 1+ failures AND time-censored data AND UPPER tail --> Get k-factor for r+1
46+
.zerofailures == FALSE & .time == TRUE & .failure == FALSE & .upper == TRUE ~ qchisq(p, df = 2*(r + 1) ) / (2*r),
47+
# 1+ failures AND time-censored data AND LOWER tail --> Get k-factor for r as normal
48+
.zerofailures == FALSE & .time == TRUE & .failure == FALSE & .upper == FALSE ~ qchisq(p, df = 2*(r) ) / (2*r),
49+
50+
# 1+ failures AND time-censored data AND LOWER tail --> Get k-factor with adjustment
51+
.zerofailures == FALSE & .time == FALSE & .failure == TRUE & .upper == TRUE ~ qchisq(p, df = 2*( (r-1) + 1)) / (2 * (r-1)) * (r-1)/r,
52+
# 1+ failures AND time-censored data AND LOWER tail --> Get k-factor with r as normal
53+
.zerofailures == FALSE & .time == FALSE & .failure == TRUE & .upper == FALSE ~ qchisq(p, df = 2*r) / (2*r),
54+
55+
# If zero failures --> then time-censored --> time = TRUE, and upper/lower distinction doesn't matter.
56+
.zerofailures == TRUE & .time == TRUE & .failure == FALSE ~ -log(1-p),
57+
# Otherwise, return NA.
58+
TRUE ~ NA_real_
59+
)
60+
61+
if(any(is.na(k))){ message("At least 1 k-factor could not be calculated, due to improper inputs. Review the rules for time-censored, failure-censored, and zero-failure data.")}
62+
63+
return(k)
64+
}
65+
66+
#' @name rk
67+
#' @title k-factor Random Deviates
68+
#' @description
69+
#' Get a random sample of k-factor values for simulating sampling distributions of failure rates.
70+
#' @param n:int number of observations.
71+
#' @param r:int number of failures (non-negative integers; can include zero)
72+
#' @param .time:logical logical; is case time-censored data?
73+
#' @param .failure:logical logical; is case failure-censored data?
74+
rk = function(n, r, .time = FALSE, .failure = FALSE){
75+
76+
# Testing values
77+
# n = 100; r = 20; .time = FALSE; .failure = FALSE
78+
79+
# Input error handling
80+
stopifnot(is.integer(as.integer(n)) & as.integer(n) > 0)
81+
stopifnot(is.logical(.time))
82+
stopifnot(is.logical(.failure))
83+
stopifnot(
84+
(.time == TRUE & .failure == FALSE) |
85+
(.time == FALSE & .failure == FALSE) |
86+
(.time == FALSE & .failure == TRUE) )
87+
88+
# Does r == 0?
89+
.zerofailures = r == 0
90+
91+
# Generate a uniform distribution of percentiles p
92+
p_uniform = runif(n = n, min = 0, max = 1)
93+
94+
# Return quantiles for the random percentiles
95+
k = qk(p = p_uniform, r = r, .time = .time, .failure = .failure)
96+
return(k)
97+
}
98+
99+
# rk(n = 1000, r = 20) %>% hist()
100+
101+
102+
#' @name pk
103+
#' @title k-factor Cumulative Distribution Function
104+
#' @description
105+
#' Function to return cumulative probabilities / percentiles given a supplied k-factor quantile `q`.
106+
#' Intended for confidence intervals for failure rates.
107+
#' @param q:[dbl] vector of quantiles (k-factors)
108+
#' @param r:int number of failures (non-negative integers; can include zero)
109+
#' @param .time:logical logical; is case time-censored data?
110+
#' @param .failure:logical logical; is case failure-censored data?
111+
pk = function(q, r, .time = FALSE, .failure = FALSE){
112+
# Testing values
113+
# q = 2; r = 20; .time = FALSE; .failure = FALSE
114+
# We just need to map the function...
115+
116+
# Construct an approximation function f, which gives the inverse of the Quantile Function,
117+
# such that you use linear interpolation to return a Probability for any Quantile supplied.
118+
by = 0.001
119+
p_range = c(by/10000, by/1000, by/100, by/10,
120+
seq(from = 0, to = 1, by = 0.001),
121+
1 - by/10, 1 - by/100, 1 - by/1000, 1 - by/10000)
122+
p_range = sort(p_range)
123+
# Get the quantiles for that range
124+
q_range = qk(p = p_range, r = r, .time = .time, .failure = .failure)
125+
# Get the inverse quantile function
126+
f = approxfun(x = q_range, y = p_range, method = "linear", rule = 2, na.rm = TRUE)
127+
# Return the expected CDF for that quantile
128+
p = f(q)
129+
return(p)
130+
}
131+
132+
133+
#' @name dk
134+
#' @title k-factor Probability Density Function
135+
#' @description
136+
#' Function to return probability densities given a supplied k-factor quantile `q`.
137+
#' Intended for visualizing sampling distributions of failure rates.
138+
#' @param x:[dbl] vector of quantiles (k-factors)
139+
#' @param r:int number of failures (non-negative integers; can include zero)
140+
#' @param .time:logical logical; is case time-censored data?
141+
#' @param .failure:logical logical; is case failure-censored data?
142+
dk = function(x, r, .time = FALSE, .failure = FALSE){
143+
# Testing values
144+
# x = 2; r = 20; .time = FALSE; .failure = FALSE
145+
146+
# Construct a range of quantiles corresponding to a range of cumulative probabilities
147+
by = 0.001
148+
p_range = c(by/10000, by/1000, by/100, by/10,
149+
seq(from = 0, to = 1, by = 0.001),
150+
1 - by/10, 1 - by/100, 1 - by/1000, 1 - by/10000)
151+
p_range = sort(p_range)
152+
# Get the quantiles for that range
153+
q_range = qk(p = p_range, r = r, .time = .time, .failure = .failure)
154+
# Fit a density curve to that quantile data, truncated at 0.
155+
curve = density(q_range, cut = c(0))
156+
157+
# Approximate a function
158+
f = approxfun(density(q_range, cut = c(0)), method = "linear", rule = 2, na.rm = TRUE)(x)
159+
# Estimate density
160+
d = f(x)
161+
return(d)
162+
}
163+
# Example
164+
# dk(x = c(0, 1, 2, 3), r = 21, .time = TRUE, .failure = FALSE)

0 commit comments

Comments
 (0)