R help - Feb 2008 - Confidence Interval for SMR

Hello,

I am looking for a function which allows to calculate the confidence
interval for a standard mortality ratio. I do have vectors with the
number of observed and expected death. Has anybody a hint where to
look?
Best,
Stefan

You could have a look at package "epitools" (on CRAN, or at
www.epitools.net).

-Heinrich.
> -----Urspr?ngliche Nachricht-----
> Von: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] Im Auftrag von stefan lhachimi
> Gesendet: Dienstag, 19. Februar 2008 10:45
> An: r-help at stat.math.ethz.ch
> Betreff: [R] Confidence Interval for SMR
> 
> 
> 
> Hello,
> 
> I am looking for a function which allows to calculate the confidence
> interval for a standard mortality ratio. I do have vectors with the
> number of observed and expected death. Has anybody a hint where to
> look?
> Best,
> Stefan
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide 
> commented, minimal, self-contained, reproducible code.
>

> Hello,
> 
> I am looking for a function which allows to calculate the confidence
> interval for a standard mortality ratio. I do have vectors with the
> number of observed and expected death. Has anybody a hint where to
> look?
> Best,
> Stefan
 We find the following function useful for an SMR
    cipoisson(observed, expected)
    
    	Terry Therneau
    	Mayo Clinic
    
"cipoisson" <- 
function(k, time = 1, p = 0.95, method = c("exact",
"anscombe") ) {
	nn <- max(length(k), length(time), length(p))
	if(nn > 1) {
		k <- rep(k, length = nn)
		time <- rep(time, length = nn)
		p <- rep(p, length = nn)
	}
	p <- (1 - p)/2  #two sided

	method <- match.arg(method)
	if(method == "exact") {
		dummy1 <- ifelse(k == 0, 1, k)
		#avoid an error message of qgamma
		lower <- ifelse(k == 0, 0, qgamma(p, dummy1))
		upper <- qgamma(1 - p, k + 1)
	}
	else if(method == "anscombe") {
		# anscombe's method
		upper <- (sqrt(k + 7/8) - qnorm(p)/2)^2
		lower <- (sqrt(k - 1/8) + qnorm(p)/2)^2
	}
	else stop("Invalid method")
	if(nn == 1)
		c(lower = lower, upper = upper)/time
	else cbind(lower = lower, upper = upper)/time
}