Abstract

For this paper, I will investigate if wearing cushioned shoes affects the incidence of ankle sprains. This project is not meant to solve a global crisis but to explore a personal interest. When not kept indoors by COVID-19, I frequently play basketball. I have twisted and sprained my ankle many times while playing the sport. Ankle sprains are not a good time. I found this dataset and corresponding study linked through the University of Florida’s stat website [References #1][References #2].

This dataset consists of 44559 observations of three variables. Each observation is an exposure to risk of injury. Exposures can be basketball practices, scrimmages, games. Naturally, the ratio of ankle sprains to exposures is low. We will only use two of these three variables in this study, SHOETYPE and SPRAIN. The variable SHOETYPE takes values of 1 and 2 to signify, respectively, if the participant is wearing cushioned or noncushioned shoes in the exposure. The variable SPRAIN takes a value of 1 if the participant sustained an ankle injury and 0 otherwise.

We summarized this data into a 2x2 contingency table. From there I used a Chi-square test to test if there if there is an association between cushioned shoes and ankle sprains. I also find the power of this test. The methods that I follow are outlined in Chapters 2 and 3 of Lachin’s Biostatistical Methods.

From these methods, we cannot reject the hypothesis that cushioned shoes have an impact on the incidence of ankle sprains. Our null hypothesis with a 95% confidence level had a p-value of .1526. The power of this test was .752, which is lower than what some experts expect. There were also some concerns that the study used convenience sampling techniques rather than random sampling techniques. All 230 basketball players were observed from 22 universities whose athletic trainers agreed to participate in the program. Considering there was still a notable decrease in ankle sprains with cushioned shoes, I would recommend modifying and repeating this experiment. If this experiment were to be repeated, I would explore the incidence of ankle sprains across a wider basketball-playing population and not just college athletes. This likely would increase the number of exposures as well, which could increase the power of the study.

The procedure and code attached to this test is as follows.

# Read-in required packages
library(tidyverse)

Procedure

First, we read in the data. We define \(n_1\), \(p_1\), \(n_2\), and \(p_2\) as the number and proportion of ankle sprains that occurred with cushioned and uncushioned shoes, respectively. We denote the pooled incidence of spraining an ankle to be \(p\) and the sample size to be N.

# Read-in data
df <- readxl::read_xlsx("ankleshoe.xlsx")
head(df)
# define p and N
p = mean(df$SPRAIN)
N = length(df$SPRAIN)

Exact 95% Asymmetric Confidence Limits for the Incidence of Ankle Sprains. We use the complementary log-log transformation of \(p\) to find these bounds. \[ \hat{\theta}=log(-log(p)) \] We then use the asymptotic normal distribution to find the 95% confidence bounds for \(\hat{\theta}\). \[ (\hat{\theta_L},\hat{\theta_U}) = \hat{\theta} \pm Z_{97.5}\sqrt{\frac{1-p}{np*log(p)^2}}\]

theta_hat <- log(-log(p))
theta_bounds = c(0,0)
theta_bounds[1] = theta_hat + qnorm(.975,0,1)*sqrt((1-p)/(N*p*log(p)^2))
theta_bounds[2] = theta_hat - qnorm(.975,0,1)*sqrt((1-p)/(N*p*log(p)^2))

Once we have these theta bounds, we can revert these back to our parameter bounds by using the inverse of the initial complementary log-log transformation. \[(\pi_L,\pi_U) = (exp[-exp(\theta_U)],exp[-exp(\theta_L)])\] Thus, our confidence bounds for \(\pi\), the incidence of ankle sprains, are:

ci_bounds = exp(-exp(theta_bounds))
print(ci_bounds)
## [1] 0.001198166 0.001926869

For this paper, we will investigate if the wearing shoes with protection affects the incidence of ankle sprains. We define \(n_1\), \(p_1\), \(n_2\), and \(p_2\) as the number and proportion of ankle sprains that occurred with cushioned and uncushioned shoes, respectively.

n1 <- length(df$SPRAIN[df$SHOETYPE==1])
p1 <- mean(df$SPRAIN[df$SHOETYPE==1])
n2 <- length(df$SPRAIN[df$SHOETYPE==2])
p2 <- mean(df$SPRAIN[df$SHOETYPE==2])

We estimate the risk difference by taking the difference of \(p_1\) and \(p_2\).

RD.est <- p2 - p1
print(RD.est)
## [1] 0.0006246895

We also find the relative risk by taking the quotient of \(p_1\) and \(p_2\).

RR.est <- p2/p1
print(RR.est)
## [1] 1.468746

Finally, we find the odds ratio of these observed proportions by taking:

OR.est <- (p2*(1-p1))/(p1*(1-p2))
print(OR.est)
## [1] 1.469665

We sum our raw data to create a 2x2 contingency table relating the amount of athletes that wore cushioned or non-cushioned shoes to the amount of athletes that experienced ankle sprains.

C <- t(table(df[,c(1,3)]))
C[c(2,1),]<-C[c(1,2),]
colnames(C) <- c("Cushioned","Noncushioned")
rownames(C) <- c("+","-")
print(C)
##       SHOETYPE
## SPRAIN Cushioned Noncushioned
##      +        41           27
##      -     30724        13767

We need to evaluate the association betweet these variables with a one-sided test. There are two options for doing this; the Fisher-Irwin exact test and the chi-squared test. The website, TowardsDataScience, details a heuristic for determining which test to use [Reference #3].

If the expected values of any cell in the contingency table is are less than 5, we should use the Fisher-Irwin exact test. Thus, we must calculate the expected values for each cell as outlined in section 2.6.2 of Lachin’s Biostatistical Methods.

p.hat <- sum(C[1,])/sum(C)
N <- sum(C)
E.11 <- sum(C[,1])*p.hat
E.12 <- sum(C[,2])*p.hat
E.21 <- sum(C[,1])*(1-p.hat)
E.22 <- sum(C[,2])*(1-p.hat)
C.EX <- C
C.EX[c(1:4)]<-c(E.11,E.21,E.12,E.22)
print(C.EX)
##       SHOETYPE
## SPRAIN   Cushioned Noncushioned
##      +    46.94944     21.05056
##      - 30718.05056  13772.94944
rm(E.11,E.12,E.21,E.22)

From this contingency table of expected values, we can see that no cell has an expected value below 5. This means we can use the Pearson’s Contingency Chi-Squared test for association. The test statistic of the Chi-Squared test is \(X^2_p=\frac{(ad-bc)^2N}{n_1n_2m_1m_2}\)

# To prevent an overflow error, this must be done in two steps. 
chi.test <- ((C[1,1]*C[2,2]-C[1,2]*C[2,1])^2*N)/(sum(C[1,])*sum(C[2,]))
chi.test <- chi.test / (sum(C[,1])*sum(C[,2])) 
chi.null <- qchisq(.950,df=1)

As \(\chi^2_0 = \;\) 3.8414588 > \(\chi^2_T = \;\) 2.4391019, we cannot reject the null hypothesis that ankle sprains have a lower incidence when cushioned shoes are worn.

We will now examine the power of this Chi-squared test with confidence level of 95%. To do this, we first find the power of the two-sided Z-test. The power of a test is defined as one minus the probability of making a type II error. As the \(\chi^2\) test is equivalent to the two-sided Z-test squared test, we can use the power for the Z-test to find the power of the \(\chi^2\) test. This approach is outlined in section 3.3.1 of Lachin’s Biostatistical Methods. The equation below outlines the process for finding the power (\(1-\beta\)) of the test.

\[Z_{1-\beta} = \frac{\sqrt{N}|\pi_1-\pi_2|-Z_{1-\alpha}\phi_0}{\phi_1};\; \phi_0 = \sqrt{N\sigma_0^2};\;\phi_1 = \sqrt{N\sigma_1^2}\] The parameters \(\pi_1\) and \(\pi_2\) are simply the proportion of sprains with uncushioned and cushions shoes, respectively. As the confidence level, \(\alpha\), is given to be 5%, all that is left to do is find \(\sigma_0\) and \(\sigma_1\), the standard deviations under the null and alternate hypotheses. The variance under the null hypothesis assumes that \(\pi_1 = \pi_2\), hence \(\sigma_0^2\) is simply: \[\sigma_0^2 = \pi(1-\pi)(\frac{1}{n_1}+\frac{1}{n_2})\] Where \(\pi\) is the incidence of sprains across the entire sample. The variance under the alternate hypothesis lends that: \[\sigma_1^2=\frac{\pi_1(1-\pi_1)}{n_1}+\frac{\pi_2(1-\pi_2)}{n_2}\] We calculate and take the square root of these variances to obtain the corresponding standard deviations in the code below:

vars <- c(0,0)
p = mean(df$SPRAIN) 
n1 = sum(C[,1])
n2 = sum(C[,2])
vars[1] = p*(1-p)*(1/n1 + 1/n2)
vars[2] = ((C[1,1]/n1)*(1-C[1,1]/n1))/n1 + ((C[1,2]/n2)*(1-C[1,2]/n2))/n2
stds  = sqrt(vars)
rm(p,vars)
print(stds)
## [1] 0.0003999901 0.0004299802

We now find \(\phi_0\) and \(\phi_1\) using these standard deviations.

phis = sqrt(N)*stds
print(phis)
## [1] 0.08443392 0.09076454

Now that all the parameters are known, we find \(Z_{1-\beta}\) and the associated power:

Z.power = (sqrt(N)*abs((C[1,1]/n1)-(C[1,2]/n2))-pnorm(.95)*phis[1])/phis[2]
print(pnorm(Z.power))
## [1] 0.7522876

The power calculated from the two-sided two proportion Z-test is equivalent to a 2x2 contingency table Chi-equare test. This power is 0.7522876. Thus, the probability of making a Type II error in this study is 0.2477124. Unfortunately, this type II error is higher than what is desired. \[ P(Type\;II\;Error) = \beta = .2477\]

Even though this experiment does not reveal that there is a significant difference between the incidence of ankle sprains between shoe type. The relative risk shows that athletes could be positively impacted by wearing cushioned shoes. It could be useful to reattempt this experiment using a larger sample size to decrease the risk of a Type II Error. To investigate this claim, we plot the power of this test against various sample sizes. We must assume that our observed sample proportions are consistent estimators of parameters \(\pi_1\) and \(\pi_2\).

Plot POWER vs N:

Z.power.fn = function(x) {
  (sqrt(x)*abs((C[1,1]/n1)-(C[1,2]/n2))-pnorm(.95)*phis[1])/phis[2]
}
tmp <- as.data.frame(curve(pnorm(Z.power.fn(x)),from=20,to=100000))
PWR.plt <- ggplot(data=tmp,aes(x=x,y=y))+
  geom_line()+
  labs(title="Power as a Function of Sample Size")+
  xlab("Sample Size")+ylab("Power")+
  ylim(0,1)+
  annotate("text", x = 20000, y = .84, label = "Desired Sample Size N \u2248 55000 if \u03B2 =.80")+
    annotate("text", x = 20000, y = .94, label = "Desired Sample Size N \u2248 89000 if \u03B2 =.90")+
  geom_hline(yintercept=.80,color="orange")+
  geom_hline(yintercept=.90,color="red")+
  theme_bw()
plot(PWR.plt)

The above chart reveals that to achieve an acceptable power of .80, we would need a sample size approximately greater than 55000. To achieve a power of .90, we would need a sample size approximately greater than 89000.

References:

#1 Dataset
http://users.stat.ufl.edu/~winner/data/ankleshoe.txt #2 Associated Study
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2386428/ #3 Fishers Test Criterion
https://towardsdatascience.com/fishers-exact-test-in-r-independence-test-for-a-small-sample-56965db48e87.