I need help With R

[Kodzo]

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

[Rick_SAS]

Why not use SAS?  I wrote a series of four articles on various aspects of computing funnel plots in SAS

I discuss the normal approximation and provide SAS code.