I need help With R

Kodzo — Mon, 08 Oct 2012 18:03:29 GMT

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

Re: I need help With R

Rick_SAS — Mon, 08 Oct 2012 20:53:17 GMT

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.

topic Re: I need help With R in SAS Programming

I need help With R

Re: I need help With R