I have the R program bellow to generate a funnel plot. I can use dummy_facid N OE ration to generate the plot. My question is how can I change the program for variance and dispersion to have the normal approximation instead of negative binomial? Can someone help??? Thanks Emmanuel ############################################################# ## Funnel Plots ## ------------------------------- ## First is the funnel.plot function ## Last is the code putting the funnel plots in a PDF funnel.plot = function(nvar, response,id=numeric(1),xlab="",ylab="",main="",pch=1) { ns = nvar ## This sets the x-axis mu = 0.04385 ## mu value doesn't really matter if we adjust for ## dispersion in variance calculations. ####### ## Calculate the dispersion of the data where ## we expect responses to be centered around 1 vars = 1/((ns*mu)^2) * ns*mu*(1-mu) dispersion = sum((response-1)^2/vars,na.rm=T)/(length(vars)) gensd = function(n) ## This function is called when we generate the confidence bands { sqrt(dispersion*(1/(n*mu))^2 * n * mu * (1-mu)) } plot(ns,response, xlab=xlab, ylab=ylab, ylim=c(-1.5,5.5),main=main,pch=pch,cex=0.8) ## Plot the data xgrid.min = min(ns,na.rm=T) xgrid.max = max(ns,na.rm=T) xgrid.spc = (xgrid.max-xgrid.min)/100 xgrid = seq(max(c(xgrid.min-xgrid.spc*10,1),na.rm=T), xgrid.max+xgrid.spc*10, xgrid.spc) ## Generate an appropriate x grid ## for confidence lines yu = 1+qnorm(0.95)*gensd(xgrid) ## 90% CI lines centered at 1 yl = 1-qnorm(0.95)*gensd(xgrid) lines(xgrid,yu,lty=5) lines(xgrid,yl,lty=5) yu = 1+qnorm(0.975)*gensd(xgrid) ## 95% CI lines centered at 1 yl = 1-qnorm(0.975)*gensd(xgrid) lines(xgrid,yu,lty=3) lines(xgrid,yl,lty=3) yu = 1+qnorm(0.995)*gensd(xgrid) ## 99% CI lines centered at 1 yl = 1-qnorm(0.995)*gensd(xgrid) lines(xgrid,yu,lty=2) lines(xgrid,yl,lty=2) abline(h=1.0,lty=2) legend("topright", inset=.005, title="Confidence Lines", c("90%","95%","99%"), lty=c(5,3,2), horiz=TRUE) ## Add legend #rect(10, -0.15, 25, 0.15, border="green",lwd=1.5,lty="dashed") # coloured zoom-in box ####################### Label Outliers #offset = (xgrid.max - xgrid.min)/45 offset = 0.0 pos.outlier = (response > (1+qnorm(0.95)*gensd(ns))) text(ns[pos.outlier]+offset,response[pos.outlier],label=id[pos.outlier],font=2,cex=0.8,adj=c(-0.3,0.2),col=rgb(.8,0,0)) neg.outlier = (response < (1-qnorm(0.95)*gensd(ns))) text(ns[neg.outlier]+offset,response[neg.outlier],label=id[neg.outlier],font=2,cex=0.8,adj=c(-0.3,0.2),col="blue") } setwd ("C:/2012 Analyses/Peds Pilot/Validation") dta <- read.table("Allmortality.csv",sep=",",header=T) attach(dta) pdf("TQIP Ped FunnelPlot.pdf",width=10,height=8) ## N_TQIP_Pt , Overall_OE (TQIP cases and O/E for Overall Mortality) funnel.plot(N,OE,dummy_facid, xlab = "Total Number of TQIP TBI Patients", ylab = "Mortality O/E Ratios", main = "Risk Adjusted Mortality: ALL Pediatric Patients") dev.off()
... View more