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()
