Untitled

Operative Technique

Sum of all data in operative technique header is 176. There are four treatments – Open, Laproscopic, lap to open, and Robotic. We want to compare the efficacy of the treatments to each other. In otherwords, we need to do 4 permutation tests where we compare a treatment to its complement. For example, we want to compare the efficacy of a laproscopic procedure to a not laproscopic procedure.

Example 1. Laproscopic vs Not Laproscopic.

Observed Statistic: In total 31 people received laproscopic procedure and 145 did not. The success rate for the laprocopic procedure, denoted as S_lap, is 22/(22+9) = 22/31 = 70.97%. The succes rate for not laproscopic, denoted S_not_lap, 90/145 = 62.07%

Our null hypothesis is that

H0: S_lap/S_not_lap = 1. H1: S_lap/S_not_lap > 1

Our observed statistic is (22/31)/(90/145).

To test, we generate a random sequence of 0s and 1s which look like our population, shuffle, and compute the above statistic. This produces a distribution. Then the p-value is the proportion of bootstraps that are equal to or less than the observed statistic.

# Step 1: Generate data

data <- c()
number_of_successes <- 99
number_of_failures <- 77

data <- c(rep(1,99), rep(0,77))

data

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

#Step 2 - Generate distribution

size_X <- 31
size_Y <- 145

propPermTest <- function(Data, num_boots, size_X, size_Y){
  #shuffle data
  output <- matrix(NA, ncol = 1, nrow = num_boots)
  p<-size_X/length(Data)
  for(i in 1:num_boots){
    #sample p% of the data
    X_index <- sample(1:length(Data), floor(p * length(Data)))
    Y_index <- setdiff(1:length(Data), X_index)
    X <- Data[X_index]
    Y <- Data[Y_index]
    #calculate the difference
    diff <- (sum(X)/size_X)/(sum(Y)/size_Y)
  
    #store
    output[i, ] <- diff
  } 
  return(output)
}

distribution <- propPermTest(data, 10000, size_X, size_Y)
length(distribution) #should be 10000

## [1] 10000

#Step 3 - plot distribution
observed_statistic <- (22/31)/(90/145)
hist(distribution)
abline(v =observed_statistic, col = "red")

#Step 4 compute p-value

sum(distribution > observed_statistic)/length(distribution)

## [1] 0.208