Amigos de R,

Eu estou envolvido num projeto que o coordenador deseja estimar SMR em diversos subgrupos. Aqui SMR (standardized mortality ratios) é definido como a quantidade de mortes observadas pela quantidade de mortes esperdas dada a probabilidade de um modelo de predição de morte. 

Eu fiz uma pequena função (abaixo) para estimar os SMRs dentre diversos subgrupos. Usei uma referencia de um inspirador para montar o ICs, mas percebi que dentre diversos trabalhos que fazem coisas parecidas a variedade de formulas para estimativas de CIs para os SMRs varia muito. Assim, como tenho muitos dados, achei que talvez fosse interessante fazer intervalos de confiança com bootstrap. Mas eu não estou sabendo ao certo como usar a função, já que a função pede pesos ou frequencias.

> dput(a[sample(1:nrow(a),100),c("Hosp_Death","SAPS3Pro2")])
structure(list(Hosp_Death = c(0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L), SAPS3Pro2 = c(0.1205, 
0.3768, 0.315, 0.3768, 0.0311, 0.0719, 0.1735, 0.75, 0.01, 0.1454, 
0.1093, 0.0511, 0.0891, 0.3353, 0.0073, 0.0401, 0.0044, 0.0239, 
0.5246, 0.1591, 0.0453, 0.0511, 0.8715, 0.1735, 0.6037, 0.4405, 
0.2047, 0.0208, 0.1591, 0.2047, 0.1205, 0.2759, 0.0273, 0.0044, 
0.2759, 0.0157, 0.0643, 0.2215, 0.2047, 0.1735, 0.1591, 0.0273, 
0.0062, 0.0311, 0.0157, 0.4192, 0.0643, 0.1735, 0.239, 0.4618, 
0.0401, 0.0802, 0.0273, 0.0239, 0.0025, 0.0157, 0.0311, 0.0157, 
0.0719, 0.239, 0.0574, 0.0157, 0.0085, 0.0085, 0.1093, 0.2215, 
0.1591, 0.0719, 0.0044, 0.8563, 0.0453, 0.2047, 0.1454, 0.0273, 
0.0891, 0.0802, 0.0239, 0.0135, 0.0802, 0.2571, 0.1735, 0.0157, 
0.3768, 0.2759, 0.0891, 0.2047, 0.1735, 0.1454, 0.0643, 0.0311, 
0.0453, 0.2759, 0.0239, 0.0181, 0.4405, 0.0085, 0.0311, 0.1205, 
0.6579, 0.1093)), .Names = c("Hosp_Death", "SAPS3Pro2"), row.names = c(21875L, 
6572L, 48643L, 11869L, 35683L, 48246L, 23919L, 29666L, 3676L, 
45549L, 36453L, 14510L, 5082L, 1898L, 41549L, 25481L, 28590L, 
38198L, 12822L, 12834L, 33267L, 34088L, 47720L, 30063L, 18326L, 
11582L, 11454L, 34960L, 18785L, 11385L, 20605L, 28105L, 25496L, 
1607L, 48766L, 36534L, 41868L, 45312L, 37206L, 26927L, 38944L, 
21087L, 22343L, 1332L, 11496L, 29485L, 38316L, 4486L, 19757L, 
45768L, 33028L, 12205L, 13150L, 41270L, 2780L, 44400L, 19696L, 
26015L, 14651L, 39093L, 24905L, 17870L, 35016L, 42851L, 20464L, 
33155L, 24924L, 33220L, 15379L, 28989L, 33286L, 34782L, 48534L, 
9045L, 20403L, 40222L, 8821L, 31240L, 12465L, 11024L, 24407L, 
45729L, 7412L, 22344L, 17737L, 3514L, 2335L, 22491L, 22493L, 
14077L, 20346L, 12114L, 36868L, 16431L, 25112L, 31272L, 10877L, 
14391L, 8422L, 45130L), class = "data.frame")

> SMR <- function(obs.var,pred.var){
+   if(length(obs.var) != length(pred.var)){
+     stop("Length of pred.var and obs.var differ.")
+   }
+   if(any(min(pred.var) <0 | max(pred.var) > 1)){
+     stop("The individual predicted death must range from 0 to 1.")
+   }
+   if(any(levels(as.factor(obs.var)) != c(0,1))){
+     stop("Observed death variable must be coded as 0 and 1.")
+   }
+   O <- length(obs.var[which(obs.var==1)])
+   E <- length(pred.var) * mean(pred.var) 
+   SMR <- O / E
+   lowerCL <- O / E * (1 - 1 / (9 * O) - 1.96 / (3 * sqrt(O)))^3 * 100
+   upperCL <- (O + 1) / E * (1 - (1 / (9 * (O + 1))) + 1.96 / (3 * sqrt(O + 1)))^3 * 100
+   output <- c(SMR=SMR,lower.Cl=lowerCL,upper.Cl=upperCL)
+   output
+ }
> SMR(a$Hosp_Death,a$SAPS3Pro2)
       SMR   lower.Cl   upper.Cl 
  1.000319  97.855941 102.244112 

O intervalo de confiança abaixo roda mas retorna valores sempre iguais a zero. 

SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}
# boot(city, ratio, R = 999, stype = "w")
boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "multicore"))
x <- boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "snow"))

Pedro Brasil