18 Set
2015
18 Set
'15
19:16
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