Differences

This shows you the differences between two versions of the page.

--- wald_test [2023/12/07 00:29] – hkimscil
+++ wald_test [2023/12/07 23:51] (current) – [Wald test in logistic regression] hkimscil
@@ Line 62: / Line 62: @@
   * Terms: A vector that specifies which coefficients to test
+====== Wald test in logistic regression ======
+<code>
+odds       <- function(p)      p/(1-p)
+odds.ratio <- function(p1, p2) odds(p1)/odds(p2)
+logit      <- function(p)      log(p/(1-p))
+ilogit     <- function(x)      exp(x)/(1+exp(x))
+iter <- 10000
+n <- 350
+p.cancer <- 0.08
+p.mutant <- 0.39
+logor <- rep (NA, iter)
+pp0 <- rep (NA, iter)
+pp1 <- rep (NA, iter)
+op0 <- rep (NA, iter)
+op1 <- rep (NA, iter)
+or <- rep (NA, iter)
+for(i in 1:iter){
+  c <- runif(n, 0, 1)
+  canc <- ifelse(c>=p.cancer, "nocancer", "cancer")
+  c <- runif(n, 0, 1)
+  gene <- ifelse(c>=p.mutant, "norm", "mutated")
+  da <- data.frame(gene, canc)
+  tab <- table(da)
+  pp0[i] <- tab[1,1] / (tab[1,1] + tab[1,2])
+  pp1[i] <- tab[2,1] / (tab[2,1] + tab[2,2])
+  op0[i] <- odds(pp0[i])
+  op1[i] <- odds(pp1[i])
+  or[i] <- odds.ratio(pp0[i], pp1[i])
+  # stats <- c(pp0, pp1, op0, op1, ortemp)
+  logor[i] <- log(or[i])
+}
+hist(logor,breaks = 50)
+</code>