Last active
August 29, 2015 14:17
-
-
Save rpetit3/0319233f6ae4100aa865 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Simulation of proposal power analysis | |
| # Y (response) --- Expect Pr(Human) >> Pr(Non-Human) | |
| # Expect N_samples > 5000 | |
| # Testing Pr(Human) = c(0.8, 0.85, 0.90, 0.95) | |
| # | |
| # X (predictor) --- Prescence or Abscence of accessory gene | |
| # Expect 2000-4000 accessory genes (tests) | |
| # Testing Pr(Accessory Gene) = c(0.05, 0.10, 0.15, 0.20, 0.25) | |
| # | |
| # Odds Ratios --- 1.1 (small), 1.5 (medium), 2.5 (large) | |
| library(foreach) | |
| library(doMC) | |
| estimate_intercept <- function(N=500, p_accessory=0.05, p_human=0.95, | |
| iterations=10) { | |
| # Estimate B0 (intercept) and B1 | |
| result <- replicate( | |
| n = iterations, | |
| expr = { | |
| X <- rbinom(N, 1, p=p_accessory) | |
| Y <- rbinom(N, 1, p=p_human) | |
| model <- glm(Y ~ X, family = "binomial") | |
| c(B0=summary(model)$coefficients[1,1], B1=summary(model)$coefficients[2,1]) | |
| } | |
| ) | |
| B0 <- sum(result[1,])/length(result[1,]) | |
| return(B0) | |
| } | |
| simulate_regression <- function(N=500, B0=0, B1=0.1, p_accessory=0.05, | |
| iterations=10) { | |
| result <- replicate( | |
| n = iterations, | |
| expr = { | |
| X1 <- rbinom(N, 1, p=p_accessory) | |
| # Fitted logit response function | |
| p_hat <- B0 + (B1 * X1) | |
| # Estimate the logit response | |
| prob <- exp(p_hat)/(1 + exp(p_hat)) | |
| uniform <- runif(length(X1), 0, 1) | |
| Y <- ifelse(uniform < prob,1,0) | |
| # Extract p-value of predictor | |
| summary(model <- glm(Y ~ X1, family = "binomial"))$coefficients[2,4] | |
| } | |
| ) | |
| # Return all p-values | |
| return(result) | |
| } | |
| # Model Parameters | |
| iterations <- 10000 | |
| N <- 6000 | |
| p_human <- c(0.8, 0.85, 0.90, 0.95) | |
| p_accessory <- c(0.05, 0.10, 0.15, 0.20, 0.25) | |
| odds_ratio <- c(1.1, 1.5, 2.5) | |
| # Use 20 cores to simulate power analysis | |
| registerDoMC(cores=20) | |
| power_analysis <- | |
| foreach(p_h=p_human, .combine='rbind') %:% | |
| foreach(p_a=p_accessory, .combine='rbind') %:% | |
| foreach(OR=odds_ratio, .combine='rbind') %dopar% { | |
| # Estimate a value for the intercept | |
| B0 <- estimate_intercept(N=N, p_accessory=p_a, | |
| p_human=p_h, iterations=iterations) | |
| # Set B1 to the odds ratio | |
| B1 <- log(OR) | |
| # Run logistic regression simulation | |
| p_vals <- simulate_regression(N=N, B0=B0, B1=B1, | |
| p_accessory=p_a, | |
| iterations=iterations) | |
| # Add values to dataframe include uncorrected power, and multiple | |
| # test corrected values (2000, 3000, 4000) | |
| data.frame( | |
| p_human=p_h, | |
| p_accessory=p_a, | |
| odds_ratio=OR, | |
| power=sum(p_vals < 0.05)/length(p_vals), | |
| power_2000=sum(p_vals < (0.05/2000))/length(p_vals), | |
| power_3000=sum(p_vals < (0.05/3000))/length(p_vals), | |
| power_4000=sum(p_vals < (0.05/4000))/length(p_vals) | |
| ) | |
| } | |
| power_analysis | |
| # N = 5000 | |
| p_human p_accessory odds_ratio power power_2000 power_3000 power_4000 | |
| 1 0.80 0.05 1.1 0.0812 0.0000 0.0000 0.0000 | |
| 2 0.80 0.05 1.5 0.6100 0.0047 0.0031 0.0023 | |
| 3 0.80 0.05 2.5 0.9984 0.4279 0.3708 0.3345 | |
| 4 0.80 0.10 1.1 0.1120 0.0001 0.0001 0.0001 | |
| 5 0.80 0.10 1.5 0.8873 0.0915 0.0724 0.0606 | |
| 6 0.80 0.10 2.5 1.0000 0.9831 0.9781 0.9738 | |
| 7 0.80 0.15 1.1 0.1525 0.0003 0.0001 0.0001 | |
| 8 0.80 0.15 1.5 0.9679 0.2647 0.2294 0.2094 | |
| 9 0.80 0.15 2.5 1.0000 0.9998 0.9994 0.9993 | |
| 10 0.80 0.20 1.1 0.1794 0.0006 0.0005 0.0003 | |
| 11 0.80 0.20 1.5 0.9910 0.4628 0.4217 0.3943 | |
| 12 0.80 0.20 2.5 1.0000 0.9999 0.9999 0.9998 | |
| 13 0.80 0.25 1.1 0.2033 0.0009 0.0005 0.0004 | |
| 14 0.80 0.25 1.5 0.9973 0.6226 0.5825 0.5518 | |
| 15 0.80 0.25 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 16 0.85 0.05 1.1 0.0691 0.0000 0.0000 0.0000 | |
| 17 0.85 0.05 1.5 0.4922 0.0006 0.0003 0.0001 | |
| 18 0.85 0.05 2.5 0.9855 0.1287 0.0948 0.0739 | |
| 19 0.85 0.10 1.1 0.0999 0.0001 0.0001 0.0001 | |
| 20 0.85 0.10 1.5 0.7993 0.0338 0.0252 0.0213 | |
| 21 0.85 0.10 2.5 0.9999 0.8647 0.8364 0.8125 | |
| 22 0.85 0.15 1.1 0.1267 0.0004 0.0002 0.0002 | |
| 23 0.85 0.15 1.5 0.9186 0.1206 0.0997 0.0862 | |
| 24 0.85 0.15 2.5 1.0000 0.9935 0.9903 0.9883 | |
| 25 0.85 0.20 1.1 0.1512 0.0002 0.0002 0.0001 | |
| 26 0.85 0.20 1.5 0.9629 0.2467 0.2153 0.1941 | |
| 27 0.85 0.20 2.5 1.0000 0.9992 0.9992 0.9992 | |
| 28 0.85 0.25 1.1 0.1661 0.0005 0.0003 0.0003 | |
| 29 0.85 0.25 1.5 0.9832 0.3819 0.3457 0.3162 | |
| 30 0.85 0.25 2.5 1.0000 0.9999 0.9999 0.9999 | |
| 31 0.90 0.05 1.1 0.0563 0.0000 0.0000 0.0000 | |
| 32 0.90 0.05 1.5 0.3432 0.0000 0.0000 0.0000 | |
| 33 0.90 0.05 2.5 0.9224 0.0000 0.0000 0.0000 | |
| 34 0.90 0.10 1.1 0.0805 0.0000 0.0000 0.0000 | |
| 35 0.90 0.10 1.5 0.6270 0.0047 0.0034 0.0028 | |
| 36 0.90 0.10 2.5 0.9983 0.4005 0.3384 0.2998 | |
| 37 0.90 0.15 1.1 0.0970 0.0001 0.0001 0.0001 | |
| 38 0.90 0.15 1.5 0.7861 0.0324 0.0255 0.0218 | |
| 39 0.90 0.15 2.5 1.0000 0.8352 0.8030 0.7747 | |
| 40 0.90 0.20 1.1 0.1216 0.0001 0.0001 0.0001 | |
| 41 0.90 0.20 1.5 0.8756 0.0760 0.0611 0.0509 | |
| 42 0.90 0.20 2.5 1.0000 0.9635 0.9536 0.9447 | |
| 43 0.90 0.25 1.1 0.1325 0.0001 0.0000 0.0000 | |
| 44 0.90 0.25 1.5 0.9263 0.1460 0.1203 0.1057 | |
| 45 0.90 0.25 2.5 1.0000 0.9931 0.9908 0.9887 | |
| 46 0.95 0.05 1.1 0.0422 0.0001 0.0000 0.0000 | |
| 47 0.95 0.05 1.5 0.1482 0.0000 0.0000 0.0000 | |
| 48 0.95 0.05 2.5 0.5674 0.0000 0.0000 0.0000 | |
| 49 0.95 0.10 1.1 0.0574 0.0000 0.0000 0.0000 | |
| 50 0.95 0.10 1.5 0.3333 0.0000 0.0000 0.0000 | |
| 51 0.95 0.10 2.5 0.9075 0.0000 0.0000 0.0000 | |
| 52 0.95 0.15 1.1 0.0723 0.0000 0.0000 0.0000 | |
| 53 0.95 0.15 1.5 0.4790 0.0006 0.0003 0.0001 | |
| 54 0.95 0.15 2.5 0.9809 0.0866 0.0561 0.0405 | |
| 55 0.95 0.20 1.1 0.0845 0.0000 0.0000 0.0000 | |
| 56 0.95 0.20 1.5 0.5896 0.0040 0.0029 0.0020 | |
| 57 0.95 0.20 2.5 0.9949 0.3316 0.2796 0.2466 | |
| 58 0.95 0.25 1.1 0.0834 0.0000 0.0000 0.0000 | |
| 59 0.95 0.25 1.5 0.6734 0.0121 0.0091 0.0074 | |
| 60 0.95 0.25 2.5 0.9985 0.5557 0.5006 0.4659 | |
| # N = 10,000 | |
| p_human p_accessory odds_ratio power power_2000 power_3000 power_4000 | |
| 1 0.80 0.05 1.1 0.1200 0.0000 0.0000 0.0000 | |
| 2 0.80 0.05 1.5 0.8990 0.0979 0.0803 0.0686 | |
| 3 0.80 0.05 2.5 1.0000 0.9872 0.9834 0.9805 | |
| 4 0.80 0.10 1.1 0.1934 0.0007 0.0004 0.0003 | |
| 5 0.80 0.10 1.5 0.9943 0.5396 0.5006 0.4747 | |
| 6 0.80 0.10 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 7 0.80 0.15 1.1 0.2624 0.0015 0.0008 0.0006 | |
| 8 0.80 0.15 1.5 0.9999 0.8478 0.8238 0.8056 | |
| 9 0.80 0.15 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 10 0.80 0.20 1.1 0.3221 0.0023 0.0019 0.0016 | |
| 11 0.80 0.20 1.5 1.0000 0.9570 0.9475 0.9407 | |
| 12 0.80 0.20 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 13 0.80 0.25 1.1 0.3667 0.0032 0.0023 0.0020 | |
| 14 0.80 0.25 1.5 1.0000 0.9883 0.9849 0.9819 | |
| 15 0.80 0.25 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 16 0.85 0.05 1.1 0.1038 0.0000 0.0000 0.0000 | |
| 17 0.85 0.05 1.5 0.8119 0.0357 0.0259 0.0217 | |
| 18 0.85 0.05 2.5 0.9999 0.8860 0.8612 0.8405 | |
| 19 0.85 0.10 1.1 0.1627 0.0003 0.0001 0.0001 | |
| 20 0.85 0.10 1.5 0.9795 0.3127 0.2788 0.2533 | |
| 21 0.85 0.10 2.5 1.0000 1.0000 0.9999 0.9998 | |
| 22 0.85 0.15 1.1 0.2192 0.0003 0.0002 0.0002 | |
| 23 0.85 0.15 1.5 0.9982 0.6386 0.5984 0.5725 | |
| 24 0.85 0.15 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 25 0.85 0.20 1.1 0.2608 0.0019 0.0010 0.0010 | |
| 26 0.85 0.20 1.5 0.9996 0.8418 0.8151 0.7968 | |
| 27 0.85 0.20 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 28 0.85 0.25 1.1 0.2956 0.0026 0.0018 0.0013 | |
| 29 0.85 0.25 1.5 0.9999 0.9302 0.9162 0.9071 | |
| 30 0.85 0.25 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 31 0.90 0.05 1.1 0.0768 0.0000 0.0000 0.0000 | |
| 32 0.90 0.05 1.5 0.6442 0.0042 0.0025 0.0020 | |
| 33 0.90 0.05 2.5 0.9981 0.4334 0.3737 0.3341 | |
| 34 0.90 0.10 1.1 0.1210 0.0003 0.0001 0.0001 | |
| 35 0.90 0.10 1.5 0.9111 0.0961 0.0805 0.0693 | |
| 36 0.90 0.10 2.5 1.0000 0.9866 0.9817 0.9776 | |
| 37 0.90 0.15 1.1 0.1632 0.0006 0.0005 0.0003 | |
| 38 0.90 0.15 1.5 0.9785 0.2991 0.2636 0.2410 | |
| 39 0.90 0.15 2.5 1.0000 0.9999 0.9997 0.9997 | |
| 40 0.90 0.20 1.1 0.1926 0.0003 0.0002 0.0002 | |
| 41 0.90 0.20 1.5 0.9933 0.5182 0.4753 0.4475 | |
| 42 0.90 0.20 2.5 1.0000 0.9999 0.9999 0.9999 | |
| 43 0.90 0.25 1.1 0.2236 0.0006 0.0005 0.0004 | |
| 44 0.90 0.25 1.5 0.9978 0.6735 0.6388 0.6120 | |
| 45 0.90 0.25 2.5 1.0000 1.0000 1.0000 1.0000 | |
| 46 0.95 0.05 1.1 0.0581 0.0000 0.0000 0.0000 | |
| 47 0.95 0.05 1.5 0.3460 0.0000 0.0000 0.0000 | |
| 48 0.95 0.05 2.5 0.9284 0.0000 0.0000 0.0000 | |
| 49 0.95 0.10 1.1 0.0750 0.0000 0.0000 0.0000 | |
| 50 0.95 0.10 1.5 0.6361 0.0041 0.0027 0.0013 | |
| 51 0.95 0.10 2.5 0.9975 0.3738 0.3149 0.2753 | |
| 52 0.95 0.15 1.1 0.1030 0.0001 0.0001 0.0001 | |
| 53 0.95 0.15 1.5 0.8032 0.0315 0.0238 0.0199 | |
| 54 0.95 0.15 2.5 1.0000 0.8428 0.8072 0.7827 | |
| 55 0.95 0.20 1.1 0.1163 0.0000 0.0000 0.0000 | |
| 56 0.95 0.20 1.5 0.8843 0.0809 0.0663 0.0561 | |
| 57 0.95 0.20 2.5 1.0000 0.9675 0.9577 0.9498 | |
| 58 0.95 0.25 1.1 0.1360 0.0001 0.0001 0.0000 | |
| 59 0.95 0.25 1.5 0.9340 0.1500 0.1249 0.1095 | |
| 60 0.95 0.25 2.5 1.0000 0.9943 0.9920 0.9905 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment