(Bayesian). A new test for multinucleoside-resistant (MNR) human immunodeficiency virus type 1 (HIV-1) variants was recently developed. The test maintains 96% sensitivity, meaning that, for those with the disease, it will correctly report “positive” for 96% of them. The test is also 98% specific, meaning that, for those without the disease, 98% will be correctly reported as “negative.” MNR HIV-1 is considered to be rare (albeit emerging), with about a .1% or .001 prevalence rate. Given the prevalence rate, sensitivity, and specificity estimates, what is the probability that an individual who is reported as positive by the new test actually has the disease? If the median cost (consider this the best point estimate) is about $100,000 per positive case total and the test itself costs $1000 per administration, what is the total first-year cost for treating 100,000 individuals
prior=0.1/100
sensitivity = 96/100
specificity = 98/100
# prior written as an odd (a:b), before it gets "updated" by Baye's factor
a = 1
b = 999
# Baye's factor calculated. P(+|with HIV)/P(+|without HIV)
Baye = sensitivity/(1-specificity)
# Multiply
a = a*Baye
# calculated proportion that P(+|with HIV)
p=(a/(a+b))
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 4.58%"
people = 100000
pos_actual = people*prior
neg_actual = people*(1-prior)
treat_cost = 100000
test_cost = 1000
true_neg = specificity*(neg_actual)
true_pos = sensitivity*(pos_actual)
false_neg = (1-specificity)*neg_actual
false_pos = (1-sensitivity)*pos_actual
cost = (test_cost*people) + ((true_pos + false_pos)*treat_cost)
print(paste0('The total cost to test and treat 100,000 people is $',cost))
## [1] "The total cost to test and treat 100,000 people is $1.1e+08"
What is the probability that, after 24 months, you received exactly 2 inspections?
p = dbinom(2,24,0.05)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 22.32%"
What is the probability that, after 24 months, you received 2 or more inspections?
p = 1-(dbinom(0,24,0.05)+dbinom(1,24,0.05))
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 33.92%"
What is the probability that your received fewer than 2 inspections?
p = dbinom(0,24,0.05)+dbinom(1,24,0.05)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 66.08%"
What is the expected number of inspections you should have received? What is the standard deviation?
# number of trials
n = 24
#probability of success in one month
p = 0.05
#mean number of inspections you should've received. aka mean
print(paste0('the expected number of inspections you should have received is ', n*p))
## [1] "the expected number of inspections you should have received is 1.2"
#standard deviation
print(paste0('standard deviation is ', sqrt(n*p*(1-p))))
## [1] "standard deviation is 1.06770782520313"
What is the probability that exactly 3 arrive in one hour?
p = dpois(3,10)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 0.76%"
What is the probability that more than 10 arrive in one hour?
# 1 - cumulative sum of probabilities from 0 to 10 arrivals in one hour.
p = 1-ppois(10,10)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 41.7%"
How many would you expect to arrive in 8 hours?
# 10 arrivals in one hour
lambda = 10
#in 8 hours we would expect 80 arrivals
print(paste0('I would expect ', 8*lambda))
## [1] "I would expect 80"
What is the standard deviation of the appropriate probability distribution?
#standard dev of poisson distribution is the root of its mean
print(paste0('standard deviation is ', sqrt(lambda)))
## [1] "standard deviation is 3.16227766016838"
If there are three family practice providers that can see 24 templated patients each day, what is the percent utilization and what are your recommendations?
The utilization rate is the ratio of the actual number of patients seen to the maximum number of patients that can be seen. The rate here is 111%. So there could be more patients by 11% than the services the family practice can provide. I would recommend increasing staffing or open up overtime opportunities to meet a reasonable utilization rate.
# 10 arrivals per hour for 8 hours is 80 patient arrivals
arrivals<-lambda*8
util_rate<-arrivals/(24*3)
print(paste0('The utilization rate is ', round(util_rate*100,2),'%'))
## [1] "The utilization rate is 111.11%"
#number of nurses sent to disney land without replacement from the 30 possible candidates if everyone was eligible equally
x = 5
# count of total nurses
m = 15
# count of total non-nurses
n = 15
# number of total people sent to disney land
k = 6
#probability of 5 nurses being sent if subordinate was innocent
p=dhyper(x,m,n,k)
print(paste0('Probability is ', round(p*100,2),'%'))
## [1] "Probability is 7.59%"
How many nurses would we have expected your subordinate to send?
#sample size
n = 6
# number of successes in pop
k = 15
# population size
N = 30
# mean of sample size of 6 should be
mu=n*(k/N)
print(paste0('I would expect ', mu))
## [1] "I would expect 3"
How many non-nurses would we have expected your subordinate to send?
I would have expected the subordinate to send the same number of nurses to non-nurses for a given sample population so for a sample of 6, 3 non-nurses would have been sent.
What is the probability that the driver will be seriously injured during the course of the year?
#use cumulative distribution function to calculate which specifies the number of trials up until success
#possible number of trials is 1200 hours in one year
trials = 1200
#probability of success of each trial
prob = 0.1/100
p = pgeom(trials,prob)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 69.93%"
In the course of 15 months?
trials = trials*(15/12)
p = pgeom(trials,prob)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 77.73%"
What is the expected number of hours that a driver will drive before being seriously injured?
# this will be the time required to guarantee a probability of an accident of 100%
print(paste0('The expected number of hours is ', (100/(prob*100)),' hours'))
## [1] "The expected number of hours is 1000 hours"
Given that a driver has driven 1200 hours, what is the probability that he or she will be injured in the next 100 hours?
p = pgeom(1300,0.001) - pgeom(1200,0.001)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 2.86%"
q=2
#probaility of failing more than 2 times with a lambda of 1 failure every 1000 hours.
p = 1-ppois(q,1)
print(paste0('probability = ', round(p*100,2),'%'))
## [1] "probability = 8.03%"
The expected value of fails is the mean, which is lambda = 1.
What is the probability that this patient will wait more than 10 minutes?
# A is probability of having waited 10 minutes
A = 1-punif(10,min=0,max=30)
print(paste0('probability = ', round((A)*100,2),'%'))
## [1] "probability = 66.67%"
If the patient has already waited 10 minutes, what is the probability that he/she will wait at least another 5 minutes prior to being seen?
# Baye's Theorem is used to "update" the prior probability given the patient has waited 10 minutes already: P(B|A) = P(A & B)/P(A). Probablilty of B (waiting another 5 min) given A (waiting 10 min) has occurred.
AandB = 1-punif(15,min=0,max=30)
B = 1-punif(5,10,30)
print(paste0('probability = ', round((AandB/A)*100,2),'%'))
## [1] "probability = 75%"
What is the expected waiting time?
min = 0
max = 30
#mean is expected wait time
print(paste0('expected wait time is ', round(((min+max)/2)),' min'))
## [1] "expected wait time is 15 min"
What is the expected failure time? What is the standard deviation?
print(paste0('expected failure time is 10'))
## [1] "expected failure time is 10"
print(paste0('standard deviation is ',round(1/10,2)))
## [1] "standard deviation is 0.1"
What is the probability that your MRI will fail after 8 years?
p = 1-pexp(8,.1)
print(paste0('probability is ',round(p,2),'%'))
## [1] "probability is 0.45%"
Now assume that you have owned the machine for 8 years. Given that you already owned the machine 8 years, what is the probability that it will fail in the next two years?
#probability A is the probability the machine has not failed in 8 years.
# Conditional probability calculated using joint probability of A and B
# divided by A already occuring.
A = 1-pexp(8,.1)
AandB = pexp(10,.1)-pexp(8,.1)
p=AandB/A
print(paste0('probability is ',round(p,2),'%'))
## [1] "probability is 0.18%"