#! /usr/bin/R sink(file="/dev/null") # to stop getting announcement text library(xtable) fisher.z <- function(r) { # function to get z-transform 1.15*log10((1+r)/(1-r)) } r.from.z <- function(z) { # inverse z-transform # get's r from z e <- exp(2*z) (e-1)/(e+1) } ci.pearson <- function(r,n,ci=.95,round=2) { # confidence interval of a Pearson correlation # defaults to 95% C.I. z <- fisher.z(r) norm <- qnorm((1-ci)/2) den <- sqrt(n-3) zl <- z + norm/den zu <- z - norm/den rl <- r.from.z(zl) ru <- r.from.z(zu) ci <- c(rl,ru) # names(ci) <- c("lower","upper") return(ci) } sink() # switch output back to stdout tag(HTML) tag(HEAD) tag(TITLE) cat("Confidence interval of observed (Pearson) correlation coefficients") untag(TITLE) untag(HEAD) lf(2) tag(BODY, bgcolor = "lime") lf(2) tag(center) cat("

Confidence interval for observed (Pearson) correlation coefficients

") comments("Let's start with some testing") prob <- 0 n <- as.numeric(scanText(formData$n)) r <- as.numeric(scanText(formData$r)) CI <- as.numeric(scanText(formData$CI)) round <- as.numeric(scanText(formData$round)) k <- length(n) k2 <- length(r) if (k != k2) { cat("

You appear to have given a different number of sample sizes and correlations, go back and try again!

") prob <- 1 } for (i in 1:k) { if ((trunc(n[i]) != n[i]) || (n[i] <= 0)) { cat("

n must be a positive integer, go back and try again!

") prob <- 1 } } for (i in 1:k) { if ((r[i] >= 1) | (r[i] < -1)) { cat("

Correlation must be between -1 and 1, go back and try again!

") prob <- 1 } } if (trunc(round) != round) { cat("

Rounding value must be an integer, go back and try again!

") prob <- 1 } if ((CI < 60) | (CI > 99.99)) { cat("

CI must be more than 60% and less than 99.99%.

") prob <- 1 } if ((round < 0) | (round > 6)) { cat("

Program only offers rounding to between zero and six decimal places. Go back and try again!

") prob <- 1 } if (prob) { cat("

Go back and try again!

") cat("CGI script written by Chris Evans using David Firth's excellent GGIwithR package. ") linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?webmaster") ; br() } lf(2) comments("end of testing") ci <- CI/100 ci.mat <- matrix(rep(NA,2*k),ncol=2) colnames(ci.mat) <- c("lwr","upr") for (i in 1:k) { ci.mat[i,] <- ci.pearson(r[i],n[i],ci,round) } res.mat <- cbind(r,n,ci.mat) #colnames(res.mat) <- c("r","n","lwr","upr") n.tot <- sum(n) z.vals <- fisher.z(r) z.mean <- mean(z.vals) norm <- abs(qnorm((1-ci)/2)) den <- sqrt(n.tot - 3) zl <- z.mean - norm/den zu <- z.mean + norm/den rl <- round(r.from.z(zl),round) ru <- round(r.from.z(zu),round) r.mean <- round(r.from.z(z.mean),round) if (!prob) { comments("Got into results section") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("") cat("
") cat("Your input") cat("
") now <- system('date +%Y-%m-%d,%T', intern=TRUE) host <- system("echo $REMOTE_ADDR", intern=TRUE) cat("Request from: ") cat("") cat(host, " at", now,"
") cat("
") cat("n values: ") cat("") cat(n) cat("
") cat("Observed correlations: ") cat("") cat(r) cat("
") cat("C.I. wanted: ") cat("") cat(CI) cat("%") cat("
") cat("Decimal places wanted: ") cat("") cat(round) cat("
") cat("Results") cat("
") m <- xtable(res.mat,caption="separate results") print(m,type="html") tag(BR) cat("
") cat("
") cat("meta-analysis") cat("
") cat("Overall r: ") cat("") cat(r.mean) cat("
") cat("CI overall for that correlation: ") cat("") cat(paste(rl,"to",ru)) cat("
") cat("Explanation of confidence intervals") cat("
A ",CI,"% CI will embrace the population proportion, i.e. the "true"") cat("proportion in the presumed infinite population from which you took your sample of ",n) cat(" in ",CI,"% of times you take a sample.") cat(" The larger the sample you take, the tighter the CI you'll see but they'll be "right" ") cat(CI,"% of the time, i.e. "wrong": not embrace the population proportion ") wrong <- 100-CI cat(wrong,"% of the time.") cat("

In significance test terms, this means that if you are comparing this with a reference value") cat(" and that reference value is not within that CI, i.e. is either smaller than ") cat(round(CI.wanted[1],round)) cat(" or else is larger than ") cat(round(CI.wanted[2],round)) cat(", you probably have a "significant" difference between your sample and that reference value.") cat("

Statistical significance is often overvalued.") cat(" All that statistical significance is telling you is that a difference ") cat("as large or larger than that was likely to happen less than one time in twenty ") cat("if you were taking samples of size ") cat(n) cat(" from an infinitely large population in which r is the reference value.") cat(" The advantage of the CI is that it tells you the precision you'd expect from random") cat(" of samples of size ") cat(n) cat(". What is generally really important in the context of the CI you have here") cat(" is how big the difference between the reference value and your observed value actually is:") cat(" is that difference clinically or managerially interesting to you or your service?

") cat("

If you are using this to compare CORE system parameters you have") cat(" with those from CORE benchmarks or descriptors, always bear in mind differences between") cat(" your or your service's setting and methods that might have contributed to any difference") cat(" you see between your proportion and that in the reference dataset; same if you don't see") cat(" a difference!") cat("

") cat("Technicalities") cat("
") cat("Confidence interval was calculated using Fisher's z transform which does assume Gaussian distributions.") cat("see Wikipaedia article http://en.wikipedia.org/wiki/Fisher_transformation") lf(2) cat("

CGI script written by Chris Evans using David Firth's excellent GGIwithR package. ") linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?webmaster") ; br() lf() cat("

") cat("Output produced at ", date()) cat("

") lf() cat("
") } untag(BODY) untag(HTML) lf() sink(paste("/home/xychris0/R/CGIwithR-log/ci-pearson.R.",as.character(now),sep="")) cat("program = ci-pearson.R ; host = ",host, "; ", now, "; n = ",n,"; r = ",r,"; CI = ",CI,"; round = ",round,"\n") sink()