[R-br] Intervalos de confiança com boot.ci

Sexta Setembro 18 16:16:37 BRT 2015

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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20150918/559d37b5/attachment.html>