Let's assume we have two groups of patients, a control group and a treatment group. They were asked a question and the answer can be yes or no. Since I come from biolgical science where I used mainly regression methods I would use logistic regression. A friend of me from the sociological sciences would analyze this data using Chi-square test or Fisher's exact test. Now, I'm wondering if there are reason to use the one or the other method.
To clarify my question I created a dataset which I analized with several methods. The R-code is below. Heare are the results (odds ratio and p-values):
odds ratio p-value comment manually (formula) 0.474 0.324 see details below chi-square - 0.524 warnings: approx. incorrect chi-square (simulated) - 0.460 fisher's exact test 0.477 0.465 logistic regression 0.474 0.324 exact what I get manually Bayes logistic regression 0.526 0.356 shrinking effect of Bayes
All results are similar. But imagine the situation if I want to publish the results performed with logistic regression and a reviewer asks me why I don't use the "classical" Chi-square test or Fisher's exact test and vice versa?
Here is how I created the data set and the analysis:
library(arm) set.seed(12345) # size of groups n <- 50 # control m <- 60 # treatment # Create data df0 <- data.frame(group = c(rep("ctrl",n), rep("treat",m)) , out = c(rbinom(n=n, 1, 0.1), rbinom(n=m,1,0.03)) ) # Tabulate with(df0,table(group, out, useNA='ifany')) out group 0 1 ctrl 42 8 treat 57 3 # Convert to matrix for using Fisher's test or Chi-square test mx0 <- with(df0,table(group, out, useNA='ifany')) # Are there expected values lower 5? round(chisq.test(mx0)$exp,1) out group 0 1 ctrl 46.4 3.6 treat 55.6 4.4 # Chi-square test chisq.test(mx0)$p.value # 0.5242444 chisq.test(mx0, simulate.p.value=TRUE)$p.value # 0.4602699 # Fisher's exact test fisher.test(mx0)$p.value # 0.4646205 fisher.test(mx0)$estimate # 0.4769107 fisher.test(mx0)$conf.int[1:2] # 0.0703001 2.6013230 # manually (OR and CI): (or <- prod(diag(mx0)) / prod(diag(apply(mx0,2,rev))) ) # 0.4736842 # log of standard error: root of (1/a+1/b+1/c+1/d): se.ln <- sqrt(sum(1/mx0)) (cil <- exp(log(or) - qnorm(0.975)*se.ln)) # 0.1074239 (ciu <- exp(log(or) + qnorm(0.975)*se.ln)) # 2.088704 (t <- abs(log(or)/se.ln)) 2*pnorm(t, lower.tail=FALSE) # 0.323628 # Logistic regression (fit1 <- summary(glm(out ~ group, data=df0, family="binomial"))$coef) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -2.1972246 0.4714045 -4.6610172 3.146504e-06 # grouptreat -0.7472144 0.7570094 -0.9870609 3.236128e-01 exp(fit1[2,"Estimate"]) # 0.4736842 (fit2 <- summary(bayesglm(out ~ group, data=df0, family="binomial"))$coef) # Estimate Std. Error z value Pr(>|z|) # (Intercept) -2.2326782 0.4645022 -4.8066039 1.535157e-06 # grouptreat -0.6415988 0.6947115 -0.9235471 3.557222e-01 exp(fit2[2,"Estimate"]) # 0.5264501