# Author: Chris Bilder # Paper: In or out? The new flagstick dilemma for putting in golf # Read in data setwd(dir = "C:\\data") set1 <- read.csv(file = "flagstick.csv") set1 # Ball speed = 1 (low), 2 (medium), 3 (high) # Entry line = 1 (center), 2 (slightly off-center), 3 (grazing) library(PropCIs) library(exact2x2) ################################################################################ # Analysis #1 # Function to calculate confidence intervals calc.ci <- function(BallSpeed, EntryLine, Rigidity = NULL, data, alpha, round.amount = 2) { # Find outcomes for factor levels # Because I use this function for both EMGA and MyGolfSpy data, I need to account for # the flagstick type that will be in the MyGolfSpy data if(!is.null(Rigidity)) { Out <- data[data$BallSpeed == BallSpeed & data$EntryLine == EntryLine & data$Rigidity == Rigidity & data$Flagstick == "Out",] In <- data[data$BallSpeed == BallSpeed & data$EntryLine == EntryLine & data$Rigidity == Rigidity & data$Flagstick == "In",] } if(is.null(Rigidity)) { Out <- data[data$BallSpeed == BallSpeed & data$EntryLine == EntryLine & data$Flagstick == "Out",] In <- data[data$BallSpeed == BallSpeed & data$EntryLine == EntryLine & data$Flagstick == "In",] Rigidity <- "Unknown" } # print(Out) # print(In) n1 <- Out$Trials n2 <- In$Trials pi.hat1 <- Out$Success/n1 pi.hat2 <- In$Success/n2 diff <- Out$Success/n1 - In$Success/n2 # Score interval - See Section 1.2 of Bilder and Loughin (2014) score.ci.diff <- diffscoreci(x1 = Out$Success, n1 = n2, x2 = In$Success, n2 = n1, conf.level = 1 - alpha) # Agresti and Caffo (Biometrics, 2000) interval AC.ci.diff <- wald2ci(x1 = Out$Success, n1 = n1, x2 = In$Success, n2 = n2, conf.level = 1 - alpha, adjust = "AC") # Fay et al. (Biometrics, 2015) interval exact.ci.diff <- binomMeld.test(x1 = Out$Success, n1 = n1, x2 = In$Success, n2 = n2, parmtype = "difference", conf.level = 1-alpha, conf.int = TRUE, alternative = "two.sided", midp = FALSE) data.frame(BallSpeed = BallSpeed, EntryLine = EntryLine, Rigidity = Rigidity, diff = diff, score.diff.low = round(score.ci.diff$conf.int[1], round.amount), score.diff.up = round(score.ci.diff$conf.int[2], round.amount), AC.diff.low = round(AC.ci.diff$conf.int[1], round.amount), AC.diff.up = round(AC.ci.diff$conf.int[2], round.amount), exact.diff.low = -round(exact.ci.diff$conf.int[2], round.amount), exact.diff.up = -round(exact.ci.diff$conf.int[1], round.amount)) } # for loops to iterate over all possible factor-level combinations save.info <- calc.ci(BallSpeed = 1, EntryLine = 1, data = set1, alpha = 0.05/9) for (i in 1:3) { for (j in 1:3) { # Divide 0.05 by 9 for Bonferroni adjustment save.it <- calc.ci(BallSpeed = i, EntryLine = j, data = set1, alpha = 0.05/9, round.amount = 2) save.info <- rbind(save.info, save.it) } } save.info[-1,] # Remove first row - used calc.ci() first as an easy way to set up structure needed for rbind() ################################################################################ # Analysis #2 library(car) library(mcprofile) # Adjust data due to 0% and 100% observed proportions const <- 0.5 Success2 <- ifelse(test = set1$Success == 0, yes = const, no = set1$Success) set1$Success2 <- ifelse(test = Success2 == 100, yes = 100 - const, no = Success2) set1 # Create qualitative factors set2 <- set1 set2$BallSpeed <- factor(set1$BallSpeed) set2$EntryLine <- factor(set1$EntryLine) set2 levels(set2$BallSpeed) levels(set2$EntryLine) # Estimate model - separate analysis showed this model had the lowest AIC mod.fit <- glm(formula = Success2/Trials ~ Flagstick * BallSpeed + EntryLine, data = set2, weights = Trials, family = binomial(link = logit)) # Warning messages regarding "non-integer" are not of concern due to the data format summary(mod.fit) Anova(mod.fit) # Calculate profile LR intervals for odds ratios OR.name <- c( "Flagstick out vs in, Ballspeed = 1", "Flagstick out vs in, Ballspeed = 2", "Flagstick out vs in, Ballspeed = 3" ) var.name <- c("int", "Flagstick", "BallSpeed2", "BallSpeed3", "EntryLine2", "EntryLine3", "Flagstick:BallSpeed2", "Flagstick:BallSpeed3") K <- matrix(data = c( 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1 ), nrow = 3, ncol = 8, byrow = TRUE, dimnames = list(OR.name, var.name)) # Estimates and intervals linear.combo <- mcprofile(object = mod.fit, CM = K) ci.log.OR <- confint(object = linear.combo, level = 0.95, adjust = "bonferroni") exp(ci.log.OR) ############### #pi_out - pi_in # pi.hat's save.pred <- predict(object = mod.fit, type = "response", se.fit = TRUE) save.pred2 <- data.frame(set2, prop = round(set2$Success2/set2$Trials,3), pi.hat = round(save.pred$fit,4), se = round(save.pred$se.fit,4)) save.pred2 # pi.hat_out - pi.hat_in calc.diff <- function(pred, EntryLine, BallSpeed) { part <- pred[pred$BallSpeed == BallSpeed & pred$EntryLine == EntryLine,] part$pi.hat[1] - part$pi.hat[2] } save.diff <- matrix(data = NA, nrow = 3, ncol = 3, dimnames = list(EntryLine = 1:3, BallSpeed = 1:3)) for(i in 1:3) { for (j in 1:3) { save.diff[i,j] <- c(calc.diff(pred = save.pred2, EntryLine = i, BallSpeed = j)) } } save.diff # Use this with results below to verify correct betas are chosen for deltaMethod() save.res <- matrix(data = NA, nrow = 3, ncol = 3, dimnames = list(EntryLine = 1:3, BallSpeed = 1:3)) parameterNames <- c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "beta7") corr.conf.level <- 1-0.05/9 # pi_out - pi_in with Entry = 1, Ball = 1 # linear.part.out11 <- "beta0 + beta1*1 + beta2*0 + beta3*0 + beta4*0 + beta5*0 + beta6*1*0 + beta7*1*0" # linear.part.in11 <- "beta0 + beta1*0 + beta2*0 + beta3*0 + beta4*0 + beta5*0 + beta6*0*0 + beta7*0*0" # g.part represents pi_out - pi_in, note that P(success) = exp(beta*x)/(1 + exp(beta*x) = 1 - (1/(1+exp(beta*x)) g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*0 + beta4*0 + beta5*0 + beta6*1*0 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*0 + beta4*0 + beta5*0 + beta6*0*0 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 1, EntryLine = 1", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[1,1] <- calc.save$Estimate # pi_out - pi_in with Entry = 2, Ball = 1 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*0 + beta4*1 + beta5*0 + beta6*1*0 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*0 + beta4*1 + beta5*0 + beta6*0*0 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 1 EntryLine = 2", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[2,1] <- calc.save$Estimate # pi_out - pi_in with Entry = 3, Ball = 1 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*0 + beta4*0 + beta5*1 + beta6*1*0 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*0 + beta4*0 + beta5*1 + beta6*0*0 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 1 EntryLine = 3", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[3,1] <- calc.save$Estimate # pi_out - pi_in with Entry = 1, Ball = 2 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*1 + beta3*0 + beta4*0 + beta5*0 + beta6*1*1 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*1 + beta3*0 + beta4*0 + beta5*0 + beta6*0*1 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 2 EntryLine = 1", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[1,2] <- calc.save$Estimate # pi_out - pi_in with Entry = 2, Ball = 2 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*1 + beta3*0 + beta4*1 + beta5*0 + beta6*1*1 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*1 + beta3*0 + beta4*1 + beta5*0 + beta6*0*1 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 2 EntryLine = 2", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[2,2] <- calc.save$Estimate # pi_out - pi_in with Entry = 3, Ball = 2 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*1 + beta3*0 + beta4*0 + beta5*1 + beta6*1*1 + beta7*1*0)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*1 + beta3*0 + beta4*0 + beta5*1 + beta6*0*1 + beta7*0*0))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 2 EntryLine = 3", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[3,2] <- calc.save$Estimate # pi_out - pi_in with Entry = 1, Ball = 3 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*1 + beta4*0 + beta5*0 + beta6*1*0 + beta7*1*1)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*1 + beta4*0 + beta5*0 + beta6*0*0 + beta7*0*1))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 3 EntryLine = 1", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[1,3] <- calc.save$Estimate # pi_out - pi_in with Entry = 2, Ball = 3 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*1 + beta4*1 + beta5*0 + beta6*0*1 + beta7*1*1)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*1 + beta4*1 + beta5*0 + beta6*0*0 + beta7*0*1))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 3 EntryLine = 2", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[2,3] <- calc.save$Estimate # pi_out - pi_in with Entry = 3, Ball = 3 g.part <- "1 - 1/(1 + exp(beta0 + beta1*1 + beta2*0 + beta3*1 + beta4*0 + beta5*1 + beta6*1*0 + beta7*1*1)) - 1 + 1/(1 + exp(beta0 + beta1*0 + beta2*0 + beta3*1 + beta4*0 + beta5*1 + beta6*0*0 + beta7*0*1))" calc.save <- deltaMethod(object = mod.fit, g = g.part, func = "BallSpeed = 3 EntryLine = 3", level = corr.conf.level, parameterNames = parameterNames) round(calc.save, 4) save.res[3,3] <- calc.save$Estimate round(save.res, 4) # Match save.diff # Using emmeans instead library(package = emmeans) calc.est1 <- emmeans(object = mod.fit, specs = ~ Flagstick + BallSpeed + EntryLine, trans = "response") confint(object = contrast(object = calc.est1, method = "revpairwise", simple = list("Flagstick"), combine = TRUE), adjust = "bonferroni", level = 0.95) ################################################################################ # MyGolfSpy analysis setwd(dir = "C:\\data") spy1 <- read.csv(file = "MyGolfSpy.csv") spy1 # Ball speed = 1 (3'), 2 (6'), 3 (9') # Entry line = 1 (center), 2 (off-center) # for loops to iterate over all possible factor-level combinations save.spy <- calc.ci(BallSpeed = 1, EntryLine = 1, Rigidity = "Standard", data = spy1, alpha = 0.05/12) for (i in 1:3) { for (j in 1:2) { for (k in c("Standard", "High")) { # print(k) # Divide 0.05 by 12 for Bonferroni adjustment save.it <- calc.ci(BallSpeed = i, EntryLine = j, Rigidity = k, data = spy1, alpha = 0.05/12, round.amount = 2) save.spy <- rbind(save.spy, save.it) } } } save.spy[-1,] # Remove first row - used calc.ci() first as an easy way to set up structure needed for rbind() #