--- title: "migraines" author: "Chris Parrish" date: "January 24, 2016" output: pdf_document --- migraines reference: - Cannon, et al., Stat2, chapter 09, example 9.6, 9.17 Import the data. {r} data <- read.csv("TMS.csv", header=TRUE) head(data) dim(data)  Odds ratio. {r} odds.painfree.TMS <- 39 / 61 odds.painfree.placebo <- 22 / 78 odds.ratio <- odds.painfree.TMS / odds.painfree.placebo odds.ratio  Logistic regresssion using glm {r} TMS.data <- as.table(matrix(c(39, 22, 61, 78), nrow=2)) dimnames(TMS.data) <- list(treatment=c("TMS", "placebo"), painfree=c("yes", "no")) TMS.data # success and failure in cols 1 and 2 Treatment <- c(1, 0) # first row is first group TMS.glm <- glm(TMS.data ~ Treatment, family=binomial) options(show.signif.stars=FALSE) summary(TMS.glm)  Alternate way of organizing the data. {r} TMS.data <- cbind(c(39, 22), c(61, 78)) dimnames(TMS.data) <- list(treatment=c("TMS", "placebo"), painfree=c("yes", "no")) TMS.data # success and failure in cols 1 and 2 Treatment <- c(1, 0) # first row is first group TMS.glm <- glm(TMS.data ~ Treatment, family=binomial) summary(TMS.glm)  LRT (= Likelihood ratio test) for utility of a simple logistic regression model $H_0 : \beta_1 = 0$ $H_a : \beta_1 \not= 0$ {r} null.deviance <- 6.8854 residual.deviance <- .7970e-14 G <- null.deviance - residual.deviance # test statistic df <- 1 - 0 p.value <- 1 - pchisq(G, df=df) p.value  CI. {r} beta1 <- coef(TMS.glm) beta1 # beta1 = log(odds ratio) exp(beta1) # exp(beta1) = odds ratio confint(TMS.glm) # CI for log(odds ratio) exp(confint(TMS.glm)) # CI for odds ratio  Probabilities of success for each group. $$P(Success \mid Treatment = TMS)$$ {r} library(boot) beta0 <- coef(TMS.glm)[1] beta1 <- coef(TMS.glm)[2] pi.TMS <- inv.logit(beta0 + beta1 * 1) pi.TMS  $$P(Success \mid Treatment = placebo)$$ {r} pi.placebo <- inv.logit(beta0 + beta1 * 0) pi.placebo