DS PS4
JL Black
2023-10-20
1. Contingency Table
| Status Quo | New Technology | |
|---|---|---|
| Test Scores Improve | 16 | 28 |
| Test Scores Same | 24 | 14 |
| Test Scores Worse | 10 | 8 |
Hypotheses
H_0 = No association tech adoption and test scores.
H_a = There is an association.
We will us TS as short hand for Test Scores and TA for Tech Adoption.
Given this H0:
First we find the expected value of each cell:
Expected Frequencies Calculation
\[ E_{ij} = \frac{\text{row total}_i \times \text{column total}_j}{\text{grand total}} \]
Applying this formula to our table, we get the following expected frequencies:
For Test Scores Improve and Status Quo: \[ E_{11} = \frac{44 \times 50}{100} = 22 \]
For Test Scores Improve and New Technology: \[ E_{12} = \frac{44 \times 50}{100} = 22 \]
For Test Scores Same and Status Quo: \[ E_{21} = \frac{38 \times 50}{100} = 19 \]
For Test Scores Same and New Technology: \[ E_{22} = \frac{38 \times 50}{100} = 19 \]
For Test Scores Worse and Status Quo: \[ E_{31} = \frac{18 \times 50}{100} = 9 \]
For Test Scores Worse and New Technology: \[ E_{32} = \frac{18 \times 50}{100} = 9 \]
Therefore, the expected frequencies table (assuming independence between technology adoption and test results) is:
| Status Quo (Expected) | New Technology (Expected) | |
|---|---|---|
| Test Scores Improve | 22 | 22 |
| Test Scores Same | 19 | 19 |
| Test Scores Worse | 9 | 9 |
Chi-Squared Test Calculation
To determine if there’s a significant association between technology adoption and test scores, we perform a chi-squared test. The test statistic is computed using the formula:
\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
Where: - \(O_{ij}\) represents the observed frequencies. - \(E_{ij}\) represents the expected frequencies under the null hypothesis of independence.
Given our observed and expected frequencies, the chi-squared test statistic is:
\[ \chi^2 = \frac{(16-22)^2}{22} + \frac{(28-22)^2}{22} + \frac{(24-19)^2}{19} + \frac{(14-19)^2}{19} + \frac{(10-9)^2}{9} + \frac{(8-9)^2}{9} = 6.1266\]
The computed \(\chi^2\) statistic follows a chi-squared distribution with degrees of freedom calculated as:
\[ \text{df} = (\text{number of rows} - 1) \times (\text{number of columns} - 1) \]
For our 3x2 table, the degrees of freedom is \((3-1) \times (2-1) = 2\).
I can now compare the computed \(\chi^2\) statistic to a chi-squared distribution with 2 degrees of freedom to determine the significance of the association.
If t > 5.991 (critical value at α=0.05 and df=2), then we reject the null hypothesis. This means there’s a significant association between technology adoption and test results at the 0.05 significance level.
If t ≤ 5.991, then we fail to reject the null hypothesis. This means that there isn’t enough evidence to suggest a significant association between technology adoption and test results at the 0.05 significance level.
Result
\[ \chi^2 = 6.1266 > 5.991 \] I reject the null hypothesis. New Tech IS correlated to test results.
2. Multiple Comparison
a. m =10
I expect to find: \(m \times \alpha = 10 \times 0.05 = 0.5\) statistically significant coefficients.
b. m = 20
I expect to find: \(m \times \alpha = 20 \times 0.05 = 1\) statistically significant coefficients.
c. m = 100
I expect to find: \(m \times \alpha = 100 \times 0.05 = 5\) statistically significant coefficients.
For a, b, and c = these are expected whether or not they are “true.” This question illustrates the importance of FWER Family Wise Error Rate.
d. Bonferroni
\((p_i ≤ α/m) = (0.0005 = 0.05/100)\) The individual p-value cutoff is \(p_i ≤ 0.0005\).
e. Benjamini-Hochberg (BH)
The Benjamini-Hochberg (BH) procedure controls the False Discovery Rate (FDR) at level \(q\):
- Order the p-values: \(P_{(1)} \leq P_{(2)} \leq ... \leq P_{(100)}\).
- Find the largest \(k\) (k is just an index to keep track of the order) such that \(P_{(k)} \leq \left(\frac{k}{100}\right)q\).
Given \(q = 0.5\):
We look for the largest \(k\) such that \(P_{(k)} \leq 0.005k\).
Assuming that the first 5 p-values (from the truly influential variables) are very close to 0, and given the distribution of the uncorrelated variables’ p-values, it’s very likely that only these 5 will pass the BH threshold.
The False Discovery Proportion (FDP) is defined as:
\[ \text{FDP} = \frac{\text{Number of false positives}}{\text{Number of variables selected}} \]
\(P_{(1)} -through- P_{(5)} = 0\) and I want to use this \(L = max {j : p_{(j)} < qj/m}\). For p 1-5 the P value 0 is less than 1/100, 2/100, 3/100, 4/100, and 5/100.
Ok - I did it in excel - now I get it. The answer is the cutoff would be at 9. For the first 9 the MAX{p-value:q*j/m} is less than 0.5. The 10’th observation goes over. Here it is in R with mechanically increasing p-values (to force them to be there Expected values not random)
# Generating mechanically increasing p-values
first_5 <- rep(0, 5)
rest_95 <- seq(0.01, 1, length.out = 95) # Generate 95 p-values from 0.01 to 1
p_values <- c(first_5, rest_95)
# Creating a dataframe for the results
Results <- data.frame(p.value = p_values)
# Rank ordering by p-values and maintaining it as a dataframe
# Note: This step is redundant since p-values are already rank-ordered, but I'm leaving it in for clarity
Results <- Results[order(Results$p.value, decreasing = FALSE), , drop = FALSE]
m <- nrow(Results)
Results$Bonferroni <- p.adjust(Results$p.value, method = "bonferroni")
Results$BH <- p.adjust(Results$p.value, method = "BH")
Results$Rank <- c(1:m)
q = 0.5 # significance threshold
library(ggplot2)
ggplot(data = Results) +
geom_col(mapping = aes(x = Rank, y = p.value)) +
geom_abline(intercept = 0, slope = q/m, linetype = "dashed") +
geom_line(mapping = aes(x = Rank, y = Bonferroni), color = "red") +
geom_line(mapping = aes(x = Rank, y = BH), color = "blue")
print(Results)## p.value Bonferroni BH Rank
## 1 0.00000000 0 0.0000000 1
## 2 0.00000000 0 0.0000000 2
## 3 0.00000000 0 0.0000000 3
## 4 0.00000000 0 0.0000000 4
## 5 0.00000000 0 0.0000000 5
## 6 0.01000000 1 0.1666667 6
## 7 0.02053191 1 0.2933131 7
## 8 0.03106383 1 0.3882979 8
## 9 0.04159574 1 0.4621749 9
## 10 0.05212766 1 0.5212766 10
## 11 0.06265957 1 0.5696325 11
## 12 0.07319149 1 0.6099291 12
## 13 0.08372340 1 0.6440262 13
## 14 0.09425532 1 0.6732523 14
## 15 0.10478723 1 0.6985816 15
## 16 0.11531915 1 0.7207447 16
## 17 0.12585106 1 0.7403004 17
## 18 0.13638298 1 0.7576832 18
## 19 0.14691489 1 0.7732363 19
## 20 0.15744681 1 0.7872340 20
## 21 0.16797872 1 0.7998987 21
## 22 0.17851064 1 0.8114120 22
## 23 0.18904255 1 0.8219241 23
## 24 0.19957447 1 0.8315603 24
## 25 0.21010638 1 0.8404255 25
## 26 0.22063830 1 0.8486088 26
## 27 0.23117021 1 0.8561860 27
## 28 0.24170213 1 0.8632219 28
## 29 0.25223404 1 0.8697726 29
## 30 0.26276596 1 0.8758865 30
## 31 0.27329787 1 0.8816060 31
## 32 0.28382979 1 0.8869681 32
## 33 0.29436170 1 0.8920052 33
## 34 0.30489362 1 0.8967459 34
## 35 0.31542553 1 0.9012158 35
## 36 0.32595745 1 0.9054374 36
## 37 0.33648936 1 0.9094307 37
## 38 0.34702128 1 0.9132139 38
## 39 0.35755319 1 0.9168031 39
## 40 0.36808511 1 0.9202128 40
## 41 0.37861702 1 0.9234561 41
## 42 0.38914894 1 0.9265451 42
## 43 0.39968085 1 0.9294904 43
## 44 0.41021277 1 0.9323017 44
## 45 0.42074468 1 0.9349882 45
## 46 0.43127660 1 0.9375578 46
## 47 0.44180851 1 0.9400181 47
## 48 0.45234043 1 0.9423759 48
## 49 0.46287234 1 0.9446374 49
## 50 0.47340426 1 0.9468085 50
## 51 0.48393617 1 0.9488945 51
## 52 0.49446809 1 0.9509002 52
## 53 0.50500000 1 0.9528302 53
## 54 0.51553191 1 0.9546887 54
## 55 0.52606383 1 0.9564797 55
## 56 0.53659574 1 0.9582067 56
## 57 0.54712766 1 0.9598731 57
## 58 0.55765957 1 0.9614820 58
## 59 0.56819149 1 0.9630364 59
## 60 0.57872340 1 0.9645390 60
## 61 0.58925532 1 0.9659923 61
## 62 0.59978723 1 0.9673988 62
## 63 0.61031915 1 0.9687606 63
## 64 0.62085106 1 0.9700798 64
## 65 0.63138298 1 0.9713584 65
## 66 0.64191489 1 0.9725983 66
## 67 0.65244681 1 0.9738012 67
## 68 0.66297872 1 0.9749687 68
## 69 0.67351064 1 0.9761024 69
## 70 0.68404255 1 0.9772036 70
## 71 0.69457447 1 0.9782739 71
## 72 0.70510638 1 0.9793144 72
## 73 0.71563830 1 0.9803264 73
## 74 0.72617021 1 0.9813111 74
## 75 0.73670213 1 0.9822695 75
## 76 0.74723404 1 0.9832027 76
## 77 0.75776596 1 0.9841116 77
## 78 0.76829787 1 0.9849973 78
## 79 0.77882979 1 0.9858605 79
## 80 0.78936170 1 0.9867021 80
## 81 0.79989362 1 0.9875230 81
## 82 0.81042553 1 0.9883238 82
## 83 0.82095745 1 0.9891054 83
## 84 0.83148936 1 0.9898683 84
## 85 0.84202128 1 0.9906133 85
## 86 0.85255319 1 0.9913409 86
## 87 0.86308511 1 0.9920518 87
## 88 0.87361702 1 0.9927466 88
## 89 0.88414894 1 0.9934258 89
## 90 0.89468085 1 0.9940898 90
## 91 0.90521277 1 0.9947393 91
## 92 0.91574468 1 0.9953747 92
## 93 0.92627660 1 0.9959963 93
## 94 0.93680851 1 0.9966048 94
## 95 0.94734043 1 0.9972004 95
## 96 0.95787234 1 0.9977837 96
## 97 0.96840426 1 0.9983549 97
## 98 0.97893617 1 0.9989145 98
## 99 0.98946809 1 0.9994627 99
## 100 1.00000000 1 1.0000000 100And here it is with actual random ordered values from 6-100. This one goes through 10 - but that’s simply due to a margin of error induced from using random draws.
# Generating the z-scores
set.seed(123) # For reproducibility
z_scores <- c(rnorm(5, -5, 0.1), rnorm(95)) # 5 variables with a z-score corresponding to a low p-value, and 95 standard normal draws
# Converting z-scores to two-tailed p-values
p_values <- 2 * (1 - pnorm(abs(z_scores)))
# Creating a dataframe for the results
Results <- data.frame(p.value = p_values)
# Rank ordering by p-values and maintaining it as a dataframe
Results <- Results[order(Results$p.value, decreasing = FALSE), , drop = FALSE]
m <- nrow(Results)
Results$Bonferroni <- p.adjust(Results$p.value, method = "bonferroni")
Results$BH <- p.adjust(Results$p.value, method = "BH")
Results$Rank <- c(1:m)
q = 0.5
ggplot(data = Results) +
geom_col(mapping = aes(x = Rank, y = p.value)) +
geom_abline(intercept = 0, slope = q/m, linetype = "dashed") +
geom_line(mapping = aes(x = Rank, y = Bonferroni), color = "red") +
geom_line(mapping = aes(x = Rank, y = BH), color = "blue")
print(Results)## p.value Bonferroni BH Rank
## 1 4.280346e-07 4.280346e-05 1.532537e-05 1
## 2 5.086584e-07 5.086584e-05 1.532537e-05 2
## 4 5.946421e-07 5.946421e-05 1.532537e-05 3
## 5 6.130147e-07 6.130147e-05 1.532537e-05 4
## 3 1.271682e-06 1.271682e-04 2.543364e-05 5
## 72 2.093421e-02 1.000000e+00 3.489035e-01 6
## 97 2.871822e-02 1.000000e+00 3.760753e-01 7
## 44 3.008603e-02 1.000000e+00 3.760753e-01 8
## 70 4.035617e-02 1.000000e+00 4.484019e-01 9
## 18 4.922736e-02 1.000000e+00 4.922736e-01 10
## 16 7.395153e-02 1.000000e+00 6.722866e-01 11
## 6 8.633331e-02 1.000000e+00 7.050949e-01 12
## 26 9.166234e-02 1.000000e+00 7.050949e-01 13
## 57 1.214412e-01 1.000000e+00 8.087525e-01 14
## 98 1.253718e-01 1.000000e+00 8.087525e-01 15
## 56 1.294004e-01 1.000000e+00 8.087525e-01 16
## 54 1.711236e-01 1.000000e+00 8.901791e-01 17
## 95 1.736235e-01 1.000000e+00 8.901791e-01 18
## 43 2.057293e-01 1.000000e+00 8.901791e-01 19
## 8 2.058494e-01 1.000000e+00 8.901791e-01 20
## 30 2.099093e-01 1.000000e+00 8.901791e-01 21
## 11 2.209214e-01 1.000000e+00 8.901791e-01 22
## 78 2.221929e-01 1.000000e+00 8.901791e-01 23
## 45 2.270619e-01 1.000000e+00 8.901791e-01 24
## 90 2.506353e-01 1.000000e+00 8.901791e-01 25
## 29 2.550633e-01 1.000000e+00 8.901791e-01 26
## 46 2.613914e-01 1.000000e+00 8.901791e-01 27
## 87 2.727118e-01 1.000000e+00 8.901791e-01 28
## 65 2.838138e-01 1.000000e+00 8.901791e-01 29
## 21 2.856000e-01 1.000000e+00 8.901791e-01 30
## 100 3.046932e-01 1.000000e+00 8.901791e-01 31
## 23 3.048895e-01 1.000000e+00 8.901791e-01 32
## 76 3.050937e-01 1.000000e+00 8.901791e-01 33
## 64 3.084046e-01 1.000000e+00 8.901791e-01 34
## 73 3.145414e-01 1.000000e+00 8.901791e-01 35
## 91 3.204645e-01 1.000000e+00 8.901791e-01 36
## 69 3.563891e-01 1.000000e+00 9.385321e-01 37
## 33 3.707199e-01 1.000000e+00 9.385321e-01 38
## 34 3.798713e-01 1.000000e+00 9.385321e-01 39
## 27 4.021503e-01 1.000000e+00 9.385321e-01 40
## 35 4.113154e-01 1.000000e+00 9.385321e-01 41
## 49 4.354114e-01 1.000000e+00 9.385321e-01 42
## 24 4.660682e-01 1.000000e+00 9.385321e-01 43
## 74 4.781999e-01 1.000000e+00 9.385321e-01 44
## 19 4.830809e-01 1.000000e+00 9.385321e-01 45
## 41 4.872390e-01 1.000000e+00 9.385321e-01 46
## 36 4.910497e-01 1.000000e+00 9.385321e-01 47
## 75 4.914474e-01 1.000000e+00 9.385321e-01 48
## 9 4.921755e-01 1.000000e+00 9.385321e-01 49
## 84 5.193313e-01 1.000000e+00 9.385321e-01 50
## 94 5.300655e-01 1.000000e+00 9.385321e-01 51
## 25 5.319453e-01 1.000000e+00 9.385321e-01 52
## 96 5.483332e-01 1.000000e+00 9.385321e-01 53
## 58 5.588075e-01 1.000000e+00 9.385321e-01 54
## 15 5.783195e-01 1.000000e+00 9.385321e-01 55
## 37 5.796352e-01 1.000000e+00 9.385321e-01 56
## 92 5.834194e-01 1.000000e+00 9.385321e-01 57
## 62 6.154400e-01 1.000000e+00 9.385321e-01 58
## 17 6.185894e-01 1.000000e+00 9.385321e-01 59
## 71 6.234044e-01 1.000000e+00 9.385321e-01 60
## 20 6.363620e-01 1.000000e+00 9.385321e-01 61
## 48 6.407465e-01 1.000000e+00 9.385321e-01 62
## 7 6.448587e-01 1.000000e+00 9.385321e-01 63
## 67 6.540018e-01 1.000000e+00 9.385321e-01 64
## 10 6.558414e-01 1.000000e+00 9.385321e-01 65
## 88 6.634307e-01 1.000000e+00 9.385321e-01 66
## 31 6.697696e-01 1.000000e+00 9.385321e-01 67
## 47 6.870329e-01 1.000000e+00 9.385321e-01 68
## 13 6.885884e-01 1.000000e+00 9.385321e-01 69
## 82 7.000297e-01 1.000000e+00 9.385321e-01 70
## 40 7.035958e-01 1.000000e+00 9.385321e-01 71
## 61 7.042130e-01 1.000000e+00 9.385321e-01 72
## 83 7.108908e-01 1.000000e+00 9.385321e-01 73
## 12 7.189864e-01 1.000000e+00 9.385321e-01 74
## 63 7.389777e-01 1.000000e+00 9.385321e-01 75
## 86 7.400539e-01 1.000000e+00 9.385321e-01 76
## 89 7.444761e-01 1.000000e+00 9.385321e-01 77
## 39 7.596331e-01 1.000000e+00 9.385321e-01 78
## 66 7.614870e-01 1.000000e+00 9.385321e-01 79
## 32 7.679393e-01 1.000000e+00 9.385321e-01 80
## 77 7.758181e-01 1.000000e+00 9.385321e-01 81
## 51 8.000221e-01 1.000000e+00 9.385321e-01 82
## 93 8.113136e-01 1.000000e+00 9.385321e-01 83
## 99 8.136652e-01 1.000000e+00 9.385321e-01 84
## 55 8.213795e-01 1.000000e+00 9.385321e-01 85
## 85 8.254922e-01 1.000000e+00 9.385321e-01 86
## 22 8.274487e-01 1.000000e+00 9.385321e-01 87
## 60 8.290333e-01 1.000000e+00 9.385321e-01 88
## 42 8.352936e-01 1.000000e+00 9.385321e-01 89
## 79 8.561294e-01 1.000000e+00 9.512549e-01 90
## 28 8.781040e-01 1.000000e+00 9.649495e-01 91
## 80 8.895360e-01 1.000000e+00 9.668870e-01 92
## 59 9.014307e-01 1.000000e+00 9.692803e-01 93
## 14 9.118680e-01 1.000000e+00 9.700723e-01 94
## 50 9.335581e-01 1.000000e+00 9.826927e-01 95
## 38 9.506331e-01 1.000000e+00 9.855151e-01 96
## 68 9.577285e-01 1.000000e+00 9.855151e-01 97
## 53 9.658048e-01 1.000000e+00 9.855151e-01 98
## 52 9.772261e-01 1.000000e+00 9.870970e-01 99
## 81 9.954009e-01 1.000000e+00 9.954009e-01 100From the first 9 rejected hypotheses, the first 5 are truly null (because their p-values are 0). So, \(V = 5\).
Given \(V = 5\) and \(R = 9\), you can calculate the FDP:
\[ FDP = \frac{V}{R} = \frac{5}{9} \]
Thus, \(FDP \approx 0.5556\) or 55.56%.
This means that approximately 55.56% of the rejected hypotheses are false positives.
3. Application
Results = read.csv("/Users/jaredblack/Documents/PhD/2023 Fall/Data Science/BHexample.csv", header = TRUE, stringsAsFactors = FALSE)
m = nrow(Results)
Results$Bonferroni = p.adjust(Results$p.value, method = "bonferroni")
Results$BH = p.adjust(Results$p.value, method = "BH")
Results$Rank = c(1:m)
q = 0.05
ggplot(data = Results) +
geom_col(mapping = aes(x = Rank, y = p.value)) +
geom_abline(slope = q/m) +
geom_line(mapping = aes(x = Rank, y = Bonferroni), color = "red") +
geom_line(mapping = aes(x = Rank, y = BH), color = "blue")
print(Results)## Dietary.variable p.value Bonferroni BH Rank
## 1 Total calories 0.001 0.025 0.0250000 1
## 2 Olive oil 0.008 0.200 0.1000000 2
## 3 Whole milk 0.039 0.975 0.2100000 3
## 4 White meat 0.041 1.000 0.2100000 4
## 5 Proteins 0.042 1.000 0.2100000 5
## 6 Nuts 0.060 1.000 0.2500000 6
## 7 Cereals and pasta 0.074 1.000 0.2642857 7
## 8 White fish 0.205 1.000 0.4910714 8
## 9 Butter 0.212 1.000 0.4910714 9
## 10 Vegetables 0.216 1.000 0.4910714 10
## 11 Skimmed milk 0.222 1.000 0.4910714 11
## 12 Red meat 0.251 1.000 0.4910714 12
## 13 Fruit 0.269 1.000 0.4910714 13
## 14 Eggs 0.275 1.000 0.4910714 14
## 15 Blue fish 0.340 1.000 0.5328125 15
## 16 Legumes 0.341 1.000 0.5328125 16
## 17 Carbohydrates 0.384 1.000 0.5647059 17
## 18 Potatoes 0.569 1.000 0.7815789 18
## 19 Bread 0.594 1.000 0.7815789 19
## 20 Fats 0.696 1.000 0.8700000 20
## 21 Sweets 0.762 1.000 0.9071429 21
## 22 Dairy products 0.940 1.000 0.9860000 22
## 23 Semi-skimmed milk 0.942 1.000 0.9860000 23
## 24 Total meat 0.975 1.000 0.9860000 24
## 25 Processed meat 0.986 1.000 0.9860000 25a.
1-7 = Olive Oil, Whole Milk, White Meat, Proteins, Nuts, and Cereals and Pasta (1 = Total Calories)
b.
the FWER = 0.05 from a. yielded a BH of 1-6, but only 1-5 actually have a p-value less than 0.05. So the FDR is 2\7 = 28.6%. If we are willing to accept up to 10% only then we would have to drop 6 and 7 leaving Olive Oil, Whole Milk, White Meat, and Protein as significant.
c.
25% FDR would allow in one more false positive (1/7 < .25) - here that’s nuts.
d.
Even if they are correlated with the outcome that does not mean they are meaningful. There could be omitted variables and the correlation could be spurious.
4. Application
a.
# Using the Bernoulli distribution to calculate the probability
# of correctly guessing the stock's movement for 10 straight days
p_success <- 0.5 # Probability of correctly guessing the movement on any single day
days <- 10 # Number of consecutive days
# Calculating the probability of success for 10 straight days
prob_correct_10_days <- p_success^days
# Displaying the result
cat("Probability of correctly guessing the stock's movement for 10 straight days if uninformed:", prob_correct_10_days)## Probability of correctly guessing the stock's movement for 10 straight days if uninformed: 0.0009765625This is well below any alpha level we might want.
b.
It depends. Pursuant to your point in c. (below) we would need to know how many other predictions she has made. If these are her only 10 - then its directionally good evidence - its enough to pursue more evidence.
c.
# Number of emails sent out
num_emails <- 10000
# Calculating the expected number of people who will receive a correct prediction sequence
expected_correct_sequences <- num_emails * prob_correct_10_days
# Displaying the result
cat("Expected number of people who will receive a correct prediction sequence:", round(expected_correct_sequences))## Expected number of people who will receive a correct prediction sequence: 10So, I could just be the unlucky one who got a false positive.
In the FDR framework, given a 0.05 threshold we can use the math we were using before and arrive at a FDP:
# Calculate FDP
FDP <- expected_correct_sequences / (expected_correct_sequences + 1)
# Displaying the result
cat("False Discovery Proportion (FDP):", round(FDP, 4))## False Discovery Proportion (FDP): 0.9071