Data 605 homework 5

1. (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?

In Step 1, we calculate the probability of testing positive (P(Positive)) using the law of total probability. In Step 2, we use Bayes’ Theorem to calculate the probability of having the disease given a positive test result (P(Disease|Positive)).

#Given values
sensitivity <- 0.96  #Sensitivity of the test (probability of a true positive)
specificity <- 0.98  #Specificity of the test (probability of a true negative)
prevalence <- 0.001  #Prevalence rate of the disease in the population

#Step 1: Calculating P(Positive) using the law of total probability
p_positive <- (sensitivity * prevalence) + ((1 - specificity) * (1 - prevalence))

#Step 2: Calculating P(Disease|Positive) using Bayes' Theorem
p_disease_given_positive <- (sensitivity * prevalence) / p_positive

#Printing the results
cat("Step 1: Probability of testing positive (P(Positive)):", p_positive, "\n")

## Step 1: Probability of testing positive (P(Positive)): 0.02094

cat("Step 2: Probability of having the disease given a positive test result (P(Disease|Positive)):", p_disease_given_positive, "\n")

## Step 2: Probability of having the disease given a positive test result (P(Disease|Positive)): 0.04584527

2. (Binomial). The probability of your organization receiving a Joint Commission inspection in any given month is .05. What is the probability that, after 24 months, you received exactly 2 inspections? What is the probability that, after 24 months, you received 2 or more inspections? What is the probability that your received fewer than 2 inspections? What is the expected number of inspections you should have received? What is the standard deviation?

For this exercise let’s use functions dbinom and pbinom to calculate probabilities for a binomial distribution.

#Given values
probability_inspection <- 0.05  #Probability of receiving an inspection in any given month
number_of_months <- 24  #Total number of months

#Step 1: Probability of receiving exactly 2 inspections after 24 months (using dbinom function)
p_exactly_2_inspections <- dbinom(2, size = number_of_months, prob = probability_inspection)

#Step 2: Probability of receiving 2 or more inspections after 24 months (using pbinom function)
p_2_or_more_inspections <- 1 - pbinom(1, size = number_of_months, prob = probability_inspection)

#Step 3: Probability of receiving fewer than 2 inspections after 24 months (using pbinom function)
p_fewer_than_2_inspections <- pbinom(1, size = number_of_months, prob = probability_inspection)

#Step 4: Expected number of inspections (mean of the binomial distribution)
expected_inspections <- number_of_months * probability_inspection

#Step 5: Standard deviation of the number of inspections (using sqrt and dbinom functions)
standard_deviation <- sqrt(number_of_months * probability_inspection * (1 - probability_inspection))

#Printing the results
cat("Step 1: Probability of exactly 2 inspections after 24 months:", p_exactly_2_inspections, "\n")

## Step 1: Probability of exactly 2 inspections after 24 months: 0.2232381

cat("Step 2: Probability of 2 or more inspections after 24 months:", p_2_or_more_inspections, "\n")

## Step 2: Probability of 2 or more inspections after 24 months: 0.3391827

cat("Step 3: Probability of fewer than 2 inspections after 24 months:", p_fewer_than_2_inspections, "\n")

## Step 3: Probability of fewer than 2 inspections after 24 months: 0.6608173

cat("Step 4: Expected number of inspections:", expected_inspections, "\n")

## Step 4: Expected number of inspections: 1.2

cat("Step 5: Standard deviation of the number of inspections:", standard_deviation, "\n")

## Step 5: Standard deviation of the number of inspections: 1.067708

3. (Poisson). You are modeling the family practice clinic and notice that patients arrive at a rate of 10 per hour. What is the probability that exactly 3 arrive in one hour? What is the probability that more than 10 arrive in one hour? How many would you expect to arrive in 8 hours? What is the standard deviation of the appropriate probability distribution? 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?

For this exercise lets use functions dpois and ppois to calculate probabilities for a Poisson distribution.

#Given values
arrival_rate_per_hour <- 10  #Patients arrive at a rate of 10 per hour
hours <- 1  #Time period is one hour
providers <- 3  #Number of family practice providers
templated_patients_per_provider <- 24  #Each provider can see 24 templated patients

#Step 1: Probability of exactly 3 patients arriving in one hour (using dpois function)
p_exactly_3_arrivals <- dpois(3, lambda = arrival_rate_per_hour * hours)

#Step 2: Probability of more than 10 patients arriving in one hour (using 1 - ppois function)
p_more_than_10_arrivals <- 1 - ppois(10, lambda = arrival_rate_per_hour * hours)

#Step 3: Expected number of arrivals in 8 hours (use lambda * time period)
expected_arrivals_8_hours <- arrival_rate_per_hour * 8

#Step 4: Standard deviation of the number of arrivals (using sqrt and dpois functions)
standard_deviation <- sqrt(arrival_rate_per_hour * hours)

#Step 5: Calculating percent utilization
percent_utilization <- (providers * templated_patients_per_provider) / expected_arrivals_8_hours * 100

#Printing the results
cat("Step 1: Probability of exactly 3 patients arriving in one hour:", p_exactly_3_arrivals, "\n")

## Step 1: Probability of exactly 3 patients arriving in one hour: 0.007566655

cat("Step 2: Probability of more than 10 patients arriving in one hour:", p_more_than_10_arrivals, "\n")

## Step 2: Probability of more than 10 patients arriving in one hour: 0.4169602

cat("Step 3: Expected number of arrivals in 8 hours:", expected_arrivals_8_hours, "\n")

## Step 3: Expected number of arrivals in 8 hours: 80

cat("Step 4: Standard deviation of the number of arrivals:", standard_deviation, "\n")

## Step 4: Standard deviation of the number of arrivals: 3.162278

cat("Step 5: Percent utilization:", percent_utilization, "%\n")

## Step 5: Percent utilization: 90 %

4. (Hypergeometric). Your subordinate with 30 supervisors was recently accused of favoring nurses. 15 of the subordinate’s workers are nurses and 15 are other than nurses. As evidence of malfeasance, the accuser stated that there were 6 company-paid trips to Disney World for which everyone was eligible. The supervisor sent 5 nurses and 1 non-nurse. If your subordinate acted innocently, what was the probability he/she would have selected five nurses for the trips? How many nurses would we have expected your subordinate to send? How many non-nurses would we have expected your subordinate to send?

For this exercise lets use the dhyper function to calculate the probability for a Hypergeometric distribution

#Given values
total_supervisors <- 30
total_nurses <- 15
total_non_nurses <- 15
total_trips <- 6
selected_nurses <- 5
selected_non_nurses <- 1

#Step 1: Probability of selecting 5 nurses out of 6 for the trips (using dhyper function)
p_selecting_5_nurses <- dhyper(selected_nurses, total_nurses, total_non_nurses, total_trips)

#Step 2: Expected number of nurses selected (using mean of hypergeometric distribution)
expected_nurses <- (total_nurses * total_trips) / total_supervisors

#Step 3: Expected number of non-nurses selected (using mean of hypergeometric distribution)
expected_non_nurses <- (total_non_nurses * total_trips) / total_supervisors

#Printing the results
cat("Step 1: Probability of selecting 5 nurses out of 6 for the trips:", p_selecting_5_nurses, "\n")

## Step 1: Probability of selecting 5 nurses out of 6 for the trips: 0.07586207

cat("Step 2: Expected number of nurses selected:", expected_nurses, "\n")

## Step 2: Expected number of nurses selected: 3

cat("Step 3: Expected number of non-nurses selected:", expected_non_nurses, "\n")

## Step 3: Expected number of non-nurses selected: 3

5. (Geometric). The probability of being seriously injured in a car crash in an unspecified location is about .1% per hour. A driver is required to traverse this area for 1200 hours in the course of a year. What is the probability that the driver will be seriously injured during the course of the year? In the course of 15 months? What is the expected number of hours that a driver will drive before being seriously injured? Given that a driver has driven 1200 hours, what is the probability that he or she will be injured in the next 100 hours?

For this exercise lets use the pgeom function to calculate probabilities for a Geometric distribution.

#Given values
probability_per_hour <- 0.001
total_hours_per_year <- 1200
hours_in_15_months <- 1200 + (15 / 12 * 24)  #Assuming an average of 24 hours per month
hours_before_injury <- 100

#Step 1: Probability of being seriously injured during the course of the year (using pgeom function)
p_injured_in_year <- pgeom(total_hours_per_year, probability_per_hour, lower.tail = TRUE)

#Step 2: Probability of being seriously injured in 15 months (using pgeom function)
p_injured_in_15_months <- pgeom(hours_in_15_months, probability_per_hour, lower.tail = TRUE)

#Step 3: Expected number of hours before being seriously injured (using mean of geometric distribution)
expected_hours_before_injury <- 1 / probability_per_hour

#Step 4: Probability of being injured in the next 100 hours given 1200 hours of driving (using pgeom function)
p_injured_in_next_100_hours <- pgeom(hours_before_injury, probability_per_hour, lower.tail = TRUE)

#Printing the results
cat("Step 1: Probability of being seriously injured during the course of the year:", p_injured_in_year, "\n")

## Step 1: Probability of being seriously injured during the course of the year: 0.6992876

cat("Step 2: Probability of being seriously injured in 15 months:", p_injured_in_15_months, "\n")

## Step 2: Probability of being seriously injured in 15 months: 0.7081794

cat("Step 3: Expected number of hours before being seriously injured:", expected_hours_before_injury, "\n")

## Step 3: Expected number of hours before being seriously injured: 1000

cat("Step 4: Probability of being injured in the next 100 hours given 1200 hours of driving:", p_injured_in_next_100_hours, "\n")

## Step 4: Probability of being injured in the next 100 hours given 1200 hours of driving: 0.09611265

6. You are working in a hospital that is running off of a primary generator which fails about once in 1000 hours. What is the probability that the generator will fail more than twice in 1000 hours? What is the expected value?

For this exercise lets use the ppois function to calculate probabilities for a Poisson distribution.

#Given values
failure_rate_per_hour <- 1 / 1000  #Probability of failure per hour
time_period <- 1000  #Total hours

#Step 1: Probability of the generator failing more than twice in 1000 hours (using ppois function)
p_failure_more_than_twice <- 1 - ppois(2, lambda = failure_rate_per_hour * time_period)

#Step 2: Expected value (mean) for the number of failures in 1000 hours (using lambda = mean of Poisson)
expected_failures <- failure_rate_per_hour * time_period

#Printing the results
cat("Step 1: Probability of the generator failing more than twice in 1000 hours:", p_failure_more_than_twice, "\n")

## Step 1: Probability of the generator failing more than twice in 1000 hours: 0.0803014

cat("Step 2: Expected value for the number of failures in 1000 hours:", expected_failures, "\n")

## Step 2: Expected value for the number of failures in 1000 hours: 1

7. A surgical patient arrives for surgery precisely at a given time. Based on previous analysis (or a lack of knowledge assumption), you know that the waiting time is uniformly distributed from 0 to 30 minutes. What is the probability that this patient will wait more than 10 minutes? 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? What is the expected waiting time?

For this exercise lets use the punif function, which calculates probabilities for a uniform distribution.

#Given values
lower_bound <- 0  #Lower bound of the uniform distribution (in minutes)
upper_bound <- 30  #Upper bound of the uniform distribution (in minutes)

#Step 1: Probability that the patient will wait more than 10 minutes
p_waiting_more_than_10_minutes <- 1 - punif(10, min = lower_bound, max = upper_bound)

#Step 2: Probability that, having already waited 10 minutes, the patient will wait at least another 5 minutes
p_waiting_at_least_another_5_minutes <- 1 - punif(15, min = lower_bound, max = upper_bound)

#Step 3: Expected waiting time (mean of the uniform distribution)
expected_waiting_time <- (lower_bound + upper_bound) / 2

#Printing the results
cat("Step 1: Probability of waiting more than 10 minutes:", p_waiting_more_than_10_minutes, "\n")

## Step 1: Probability of waiting more than 10 minutes: 0.6666667

cat("Step 2: Probability of waiting at least another 5 minutes after already waiting 10 minutes:", p_waiting_at_least_another_5_minutes, "\n")

## Step 2: Probability of waiting at least another 5 minutes after already waiting 10 minutes: 0.5

cat("Step 3: Expected waiting time:", expected_waiting_time, "minutes\n")

## Step 3: Expected waiting time: 15 minutes

8. Your hospital owns an old MRI, which has a manufacturer’s lifetime of about 10 years (expected value). Based on previous studies, we know that the failure of most MRIs obeys an exponential distribution. What is the expected failure time? What is the standard deviation? What is the probability that your MRI will fail after 8 years? 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?

For this exercise lets use the pexp function, which calculates probabilities for an exponential distribution.

#Given values
expected_lifetime <- 10  #Expected lifetime of the MRI in years

#Step 1: Expected failure time (mean of the exponential distribution)
expected_failure_time <- expected_lifetime

#Step 2: Standard deviation of the exponential distribution
standard_deviation <- 1  #For an exponential distribution, the standard deviation is equal to the mean

#Step 3: Probability that the MRI will fail after 8 years
p_failure_after_8_years <- pexp(8, rate = 1/expected_lifetime)

#Step 4: Given that you already owned the machine for 8 years, probability of failure in the next two years
p_failure_in_next_two_years <- pexp(2, rate = 1/expected_lifetime)

#Printing the results
cat("Step 1: Expected failure time:", expected_failure_time, "years\n")

## Step 1: Expected failure time: 10 years

cat("Step 2: Standard deviation of the exponential distribution:", standard_deviation, "years\n")

## Step 2: Standard deviation of the exponential distribution: 1 years

cat("Step 3: Probability of failure after 8 years:", p_failure_after_8_years, "\n")

## Step 3: Probability of failure after 8 years: 0.550671

cat("Step 4: Probability of failure in the next two years, given ownership for 8 years:", p_failure_in_next_two_years, "\n")

## Step 4: Probability of failure in the next two years, given ownership for 8 years: 0.1812692