#modifications by John Christie #03/10/2011 # The variable "slices" in each function is there so that one can adjust the quality of the measure of the area under the curve. It's one of the main reasons that this code will give different results from online calculators. It can also be lowered in order to overcome any possible performance issues (unwise in the cross_2â€¦ functions). Everything should be fast enough as is. Bf<-function(sd, obtained, uniform = FALSE, lower=0, upper=1, meanoftheory=0, sdtheory=1, tails=2) { #Version 2.0.1 #modifications by John Christie # modification means that this does not exactly replicate Dienes. What it does # is do the right thing instead. :) The current version is more accurate. #03/10/2011 #Authors Danny Kaye & Thom Baguley #Version 1.0 #19/10/2009 # test data can be found starting at p100 # notes on how to use these functions can be found at # http://danny-kaye.co.uk/Docs/Dienes_notes.pdf #raised from 2000 for better accuracy - the speed of the new code allows it slices <- 20000 if(uniform){ range <- upper - lower dist_theta <- 1 / range # incr <- range / slices # theta <- seq(lower + incr, by = incr, length.out = slices+1) # height <- dist_theta * dnorm(obtained, theta, sd) # area <- sum(height * incr) # the commented code above replicates the original result but at the # limit (slices <- 5e6) it's actually equivalent to the followingâ€¦ area <- dist_theta * diff(pnorm(c(lower, upper), obtained, sd)) }else{ # the code below again doesnt' replicate the original code. # incrementing in a scalar loop was causing an accumulation of tiny fp errors # the lower end was incremented prior to starting (that's fixed above too) zlim <- 5 incr <- sdtheory / (slices/(zlim*2)) newLower <- meanoftheory - zlim * sdtheory theta <- seq(newLower, by = incr, length.out = slices+1) dist_theta <- dnorm(theta, meanoftheory, sdtheory) if (tails == 1){ dist_theta <- dist_theta[theta > 0] * 2 theta <- theta[theta > 0] } height <- dist_theta * dnorm(obtained, theta, sd) area <- sum(height * incr) } LikelihoodTheory <- area Likelihoodnull <- dnorm(obtained, 0, sd) BayesFactor <- LikelihoodTheory / Likelihoodnull return( list("LikelihoodTheory" = LikelihoodTheory, "Likelihoodnull" = Likelihoodnull, "BayesFactor" = BayesFactor) ) } contrast<-function(SEd, meand, nu, theta1, theta2) { #modifications by John Christie #01/10/2011 #Authors Danny Kaye & Thom Baguley #19/10/2009 #The test data is found in Box 3.7 pp 74-5 #SEd - contrast standard error #meand - contrast mean #nu - contrast degrees of freedom # slices <- 10000 zlim <- 5 inc <- SEd / (slices/(zlim*2)) theta <- seq(meand - zlim * SEd, by = inc, length.out = slices+1) likelihood <- (1 + (meand - theta)^2 / (nu * SEd^2))^(-(nu + 1) / 2) likelihoodmax <- max(likelihood) likelihood <- likelihood / likelihoodmax which32 <- which(likelihood >= 1/32) which8 <- which(likelihood >= 1/8) begin32 <- theta[which32[1]] begin8 <- theta[which8[1]] end8 <- theta[tail(which8,1)+1] end32 <- theta[tail(which32,1)+1] B1 <- round((theta1 - theta[1]) / inc + 1, 0) B2 <- round((theta2 - theta[1]) / inc + 1, 0) likelihoodratio <- likelihood[B1] / likelihood[B2] return( list("begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) ) } meanUV<-function(SEd, meand, n, theta1, theta2) { #modifications by John Christie #01/10/2011 #Authors Danny Kaye & Thom Baguley #Version 1.0 #19/10/2009 #The test data is from page 138 #SEd - sample standard error #meand - sample mean #n - number of subjects #theta1 - population mean assumed by first hypothesis #theta2 - population mean assumed by second hypothesis # Vard <- n * SEd^2 SSd <- Vard * (n - 1) slices <- 10000 zlim <- 5 inc <- SEd / (slices/(zlim*2)) theta <- seq(meand - zlim * SEd, by = inc, length.out = slices+1) likelihood <- (SSd + n * (meand - theta)^2)^(-(n - 2) / 2) likelihoodmax <- max(likelihood) likelihood <-likelihood / likelihoodmax #unlike the original I find a chunk that's the middle of the distribution which32 <- which(likelihood >= 1/32) which8 <- which(likelihood >= 1/8) begin32 <- theta[which32[1]] begin8 <- theta[which8[1]] end8 <- theta[ tail(which8, 1) + 1 ] end32 <- theta[ tail(which32, 1) + 1 ] B1 <- round((theta1 - theta[1]) / inc + 1, 0) B2 <- round((theta2 - theta[1]) / inc + 1, 0) likelihoodratio <- likelihood[B1] / likelihood[B2] return( list("begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) ) } cross2x2_odds <- function(a, b, c, d, psi1, psi2) { #modifications by John Christie #01/10/2011 #Authors Danny Kaye & Thom Baguley #Version 1.0 #19/10/2009 #The test data is from figure 5.4 p 135 # m <- a + b n <- c + d k <- a + c l <- b + d slices <- 10000 psi <- (1:slices) / sqrt(slices) starting <- 0 if((a - d) > 0) starting <- a - d ending <- m if((a + c) < m) ending <- a + c A <- starting:ending likelihood <- sapply( 1:slices, function(B) { sum( choose(m, A) * choose(n, (a + c - A)) * psi[B]^(A - a))}) likelihood <- likelihood^-1 likelihoodmax <- max(likelihood) psimax <- psi[which.max(likelihood)] likelihood <- likelihood / likelihoodmax which32 <- which(likelihood >= 1/32) which8 <- which(likelihood >= 1/8) begin32 <- psi[which32[1]] begin8 <- psi[which8[1]] end8 <- psi[ tail(which8, 1) + 1 ] end32 <- psi[ tail(which32, 1) + 1 ] psi1 <- psi1 * 100 psi2 <- psi2 * 100 likelihoodratio <- likelihood[psi1] / likelihood[psi2] return( list("psimax"=psimax, "begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) ) } cross2x22<-function(a, b, c, d, gamma1, gamma2) { #modifications by John Christie #01/10/2011 #Authors Danny Kaye & Thom Baguley #Version 1.0 #19/10/2009 #The test data is from figure 5.4 p 135 # m <- a + b n <- c + d k <- a + c l <- b + d slices <- 10000 gammamax <- 0 gamma <- (1:slices) / sqrt(slices) A <- b:(m + n - d) likelihood <- sapply( 1:slices, function(B){ sum( choose(A - 1, b - 1) * choose(m + n - A - 1, d - 1) * gamma[B]^(A - m) ) } ) likelihood <- likelihood^-1 likelihoodmax <- max(likelihood) gammamax <- gamma[which.max(likelihood)] likelihood <- likelihood / likelihoodmax which32 <- which(likelihood >= 1/32) which8 <- which(likelihood >= 1/8) begin32 <- gamma[which32[1]] begin8 <- gamma[which8[1]] end8 <- gamma[ tail(which8, 1) + 1 ] end32 <- gamma[ tail(which32, 1) + 1 ] gamma1 <- gamma1 * 10 gamma2 <- gamma2 * 10 likelihoodratio <- likelihood[gamma1] / likelihood[gamma2] return( list("gammamax"=gammamax, "begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) ) } proportion <- function(suc, fail) { #modifications by John Christie #01/10/2011 #Authors Danny Kaye & Thom Baguley #Version 1.0 #19/10/2009 #The test data is from figure 5.3 p131 #not sure this does what he says it does, no theta input is asked for #in the book code, as implied in the text and no ratio is returned #note that the function is called with successes and failures not #successes and trials # slices <- 10000 theta <- (1:slices)/slices likelihood <- theta^suc*(1-theta)^fail likelihoodmax <- max(likelihood) thetamax <- theta[which.max(likelihood)] likelihood <- likelihood / likelihoodmax #unlike the original I find a chunk that's the middle of the distribution which32 <- which(likelihood >= 1/32) which8 <- which(likelihood >= 1/8) begin32 <- theta[which32[1]] begin8 <- theta[which8[1]] end8 <- theta[ tail(which8, 1) + 1 ] end32 <- theta[ tail(which32, 1) + 1 ] return( list("thetamax"=thetamax, "begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32) ) }