Intervalos de confiança com boot.ci
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
Tentei rodar e não consegui.
dput(a[sample(1:nrow(a),100),c("Hosp_Death","SAPS3Pro2")])
Error in dput(a[sample(1:nrow(a), 100), c("Hosp_Death", "SAPS3Pro2")]) : object 'a' not found 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 --- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
vc pode fazer assim: SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2) SMR2(a) ### SMR observada nos dados require(boot) plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
Prezados, Também tentei rodar a rotina e o objeto a não foi encontrado, sedo o erro Erro em NROW(data) : objeto 'a' não encontrado -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ====================================================================== Em 22/09/2015 05:19, Elias Teixeira Krainski escreveu:
vc pode fazer assim:
SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
O objeto a é o banco de dados que precisa ser criado a partir da saida do dput(). Pedro Brasil Em 22 de setembro de 2015 11:19, ASANTOS <alexandresantosbr@yahoo.com.br> escreveu:
Prezados,
Também tentei rodar a rotina e o objeto a não foi encontrado, sedo o erro
Erro em NROW(data) : objeto 'a' não encontrado
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Em 22/09/2015 05:19, Elias Teixeira Krainski escreveu:
vc pode fazer assim:
SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
Alguém poderia passar o código anterior, pois eu apaguei o que chegou antes. vc pode fazer assim: SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2) SMR2(a) ### SMR observada nos dados require(boot) plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo) . --- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
Mauro, segue o script, 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 } SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w} boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "multicore")) x <- boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "snow")) SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2) SMR2(a) ### SMR observada nos dados require(boot) plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo) Em 22/09/2015 23:14, Mauro Sznelwar escreveu:
Alguém poderia passar o código anterior, pois eu apaguei o que chegou antes.
vc pode fazer assim:
SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
.
--- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Amigos de R, Desculpa pela demora, so consegui voltar nessa análise hoje. Bom tentei fazer o script recomendado pelo Elias. O pedaço do plot do bootstrap funciona, mas o pedaço para gerar o intervalo de confiança retorna um erro.
library(boot)> SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}> SMR2(a)[1] 1.000319> plot(boo <- boot(a,SMR2,999,stype='i'))> boot.ci(boo)Error in bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o, : estimated adjustment 'w' is infinite
Achei muito estranho a função plot do bot funcionar e o boot.ci não funcionar. Possivelmente ha necessidade de resolver alguma configuração para esse pedação do script. Pedro Brasil 2015-09-23 0:19 GMT-03:00 ASANTOS <alexandresantosbr@yahoo.com.br>:
Mauro, segue o script,
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 }
SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}
boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "multicore"))
x <- boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "snow")) SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
Em 22/09/2015 23:14, Mauro Sznelwar escreveu:
Alguém poderia passar o código anterior, pois eu apaguei o que chegou antes.
vc pode fazer assim:
SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
.
--- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
Considerando a funcao que calcula a estatistica, e' mais apropriado usar
stype='w' na funcao boot()
O CI fica menos amplo que o obtido por GLM:
plot(boo <- boot(a,SMR2,999,stype='w'))
boot.ci(boo)
m <- glm(Hosp_Death ~ offset(log(SAPS3Pro2)), poisson, data=a)
(coef.stats <- coef(summary(m)))
hist(boo$t)
abline(v=exp(coef.stats[1,1] +
qnorm(c(0.025,0.5,0.975))*coef.stats[1,2]), col=2)
abline(v=quantile(boo$t, c(0.025,0.975)), col=4)
Elias
On 23/09/15 15:40, Pedro Emmanuel Alvarenga Americano do Brasil wrote:
> Amigos de R,
>
> Desculpa pela demora, so consegui voltar nessa análise hoje. Bom
> tentei fazer o script recomendado pelo Elias. O pedaço do plot do
> bootstrap funciona, mas o pedaço para gerar o intervalo de confiança
> retorna um erro.
>
> >library(boot) > SMR2 <-
> function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2)
> * mean(a$SAPS3Pro2))*w} > SMR2(a) [1] 1.000319
> >plot(boo <- boot(a,SMR2,999,stype='i')) > boot.ci
> <http://boot.ci>(boo) Error in bca.ci <http://bca.ci>(boot.out, conf,
> index[1L], L = L, t = t.o, t0 = t0.o, : estimated adjustment 'w' is
> infinite
>
> Achei muito estranho a função plot do bot funcionar e o boot.ci
> <http://boot.ci> não funcionar. Possivelmente ha necessidade de
> resolver alguma configuração para esse pedação do script.
>
> Pedro Brasil
>
> 2015-09-23 0:19 GMT-03:00 ASANTOS <alexandresantosbr@yahoo.com.br
> <mailto:alexandresantosbr@yahoo.com.br>>:
>
> Mauro, segue o script,
>
>
> 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
> }
>
>
> SMR2 <-
> function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2)
> * mean(a$SAPS3Pro2))*w}
>
>
> boot.ci <http://boot.ci>(boot(a, SMR2, R = 10, stype =
> "w",parallel = "multicore"))
>
> x <- boot.ci <http://boot.ci>(boot(a, SMR2, R = 10, stype =
> "w",parallel = "snow"))
> SMR2 <- function(a,w=1)
> sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
>
> SMR2(a) ### SMR observada nos dados
>
> require(boot)
>
> plot(boo <- boot(a,SMR2,999,stype='i'))
> boot.ci <http://boot.ci>(boo)
>
> Em 22/09/2015 23:14, Mauro Sznelwar escreveu:
>
> Alguém poderia passar o código anterior, pois eu apaguei o que
> chegou antes.
>
>
>
> vc pode fazer assim:
>
> SMR2 <- function(a,w=1)
> sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
>
> SMR2(a) ### SMR observada nos dados
>
> require(boot)
>
> plot(boo <- boot(a,SMR2,999,stype='i'))
> boot.ci <http://boot.ci>(boo)
>
> .
>
>
> ---
> Este email foi escaneado pelo Avast antivírus.
> https://www.avast.com/antivirus
>
> _______________________________________________
> R-br mailing list
> R-br@listas.c3sl.ufpr.br <mailto:R-br@listas.c3sl.ufpr.br>
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e
> fornea cdigo mnimo reproduzvel.
>
>
> --
> ======================================================================
> Alexandre dos Santos
> Proteção Florestal
> IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato
> Grosso
> Campus Cáceres
> Caixa Postal 244
> Avenida dos Ramires, s/n
> Bairro: Distrito Industrial
> Cáceres - MT CEP: 78.200-000
> Fone: (+55) 65 8132-8112 <tel:%28%2B55%29%2065%208132-8112> (TIM)
> (+55) 65 9686-6970 <tel:%28%2B55%29%2065%209686-6970> (VIVO)
> e-mails:alexandresantosbr@yahoo.com.br
> <mailto:e-mails%3Aalexandresantosbr@yahoo.com.br>
> alexandre.santos@cas.ifmt.edu.br
> <mailto:alexandre.santos@cas.ifmt.edu.br>
> Lattes: http://lattes.cnpq.br/1360403201088680
> ======================================================================
>
> _______________________________________________
> R-br mailing list
> R-br@listas.c3sl.ufpr.br <mailto:R-br@listas.c3sl.ufpr.br>
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e
> fornea cdigo mnimo reproduzvel.
>
>
>
>
> _______________________________________________
> R-br mailing list
> R-br@listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
Elias, Na verdade o w entrou na função apenas porque eu estava tentando reproduzir o primeiro exemplo da função boot.ci. Eu não tinha intenção de colocar qualquer peso na função. Mas sem o w a função boot retonrna sempre um erro. Assim coloquei o w sempre com peso 1, achando que assim daria certo mas sempre com pesos sem "função". O fato é que agora eu voltei a ter o problema inicial. A função funciona mas retorna intervalos com valores iguais a zero. Eu esperava que o intervalo contivesse o valor que retorna a função originalmente que é 1.000319 e que fosse alguma coisa proxima dos intervalos que retornam a função original SMR. Assim, os intervalos ainda não fazem sentido para mim. Eu ainda não consegui entender o que estou fazendo de errado.
SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}> SMR2(a)[1] 1.000319> plot(boo <- boot(a,SMR2,999,stype='w'))> boot.ci(boo)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates
CALL : boot.ci(boot.out = boo) Intervals : Level Normal Basic Studentized 95% ( 0, 0 ) ( 0, 0 ) ( 0, 0 ) Level Percentile BCa 95% ( 0, 0 ) ( 0, 0 ) Calculations and Intervals on Original Scale> quantile(boo$t, c(0.025,0.975)) 2.5% 97.5% 0.000000e+00 3.777955e-09
Pedro Brasil Em 23 de setembro de 2015 10:50, Elias Teixeira Krainski < eliaskrainski@yahoo.com.br> escreveu:
Considerando a funcao que calcula a estatistica, e' mais apropriado usar stype='w' na funcao boot()
O CI fica menos amplo que o obtido por GLM:
plot(boo <- boot(a,SMR2,999,stype='w')) boot.ci(boo)
m <- glm(Hosp_Death ~ offset(log(SAPS3Pro2)), poisson, data=a) (coef.stats <- coef(summary(m)))
hist(boo$t) abline(v=exp(coef.stats[1,1] + qnorm(c(0.025,0.5,0.975))*coef.stats[1,2]), col=2) abline(v=quantile(boo$t, c(0.025,0.975)), col=4)
Elias
On 23/09/15 15:40, Pedro Emmanuel Alvarenga Americano do Brasil wrote:
Amigos de R,
Desculpa pela demora, so consegui voltar nessa análise hoje. Bom tentei fazer o script recomendado pelo Elias. O pedaço do plot do bootstrap funciona, mas o pedaço para gerar o intervalo de confiança retorna um erro.
library(boot)> SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}> SMR2(a)[1] 1.000319> plot(boo <- boot(a,SMR2,999,stype='i'))> boot.ci(boo)Error in bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o, : estimated adjustment 'w' is infinite
Achei muito estranho a função plot do bot funcionar e o boot.ci não funcionar. Possivelmente ha necessidade de resolver alguma configuração para esse pedação do script.
Pedro Brasil
2015-09-23 0:19 GMT-03:00 ASANTOS <alexandresantosbr@yahoo.com.br>:
Mauro, segue o script,
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 }
SMR2 <- function(a,w=1){length(a$Hosp_Death[which(a$Hosp_Death==1)])*w/(length(a$SAPS3Pro2) * mean(a$SAPS3Pro2))*w}
boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "multicore"))
x <- boot.ci(boot(a, SMR2, R = 10, stype = "w",parallel = "snow")) SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
Em 22/09/2015 23:14, Mauro Sznelwar escreveu:
Alguém poderia passar o código anterior, pois eu apaguei o que chegou antes.
vc pode fazer assim:
SMR2 <- function(a,w=1) sum(w*a$Hosp_Death)/sum(w*a$SAPS3Pro2)
SMR2(a) ### SMR observada nos dados
require(boot)
plot(boo <- boot(a,SMR2,999,stype='i')) boot.ci(boo)
.
--- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo mnimo reproduzvel.
_______________________________________________ R-br mailing listR-br@listas.c3sl.ufpr.brhttps://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
O argumento 'w' e' mudado internamente na funcao boot(), ou seja, nao e' usado como sendo igual a 1 internamente. Essa alteracao e' o truque usado para o processo de reamostragem considerando estratos. Creio que vc nao precisa usar o pacote boot() para fazer reamostragem. E pode fazer reamostragem simples, pois nao vi situacao de haver estratos nos seus dados. Assim, voce poderia fazer smr.fun <- function(x) sum(x$Hosp_Death)/sum(x$SAPS3Pro2) smr0 <- smr.fun(a) smr.boot <- sapply(1:9999, function(x) smr.fun(a[sample(1:nrow(a), replace=TRUE),])) ### HPD interval (Intervalo de mais alta densidade) require(coda) hpdic <- HPDinterval(as.mcmc(smr.boot), 0.95) hist(smr.boot, breaks=pretty(smr.boot, 50)) abline(v=hpdic) Elias
Elias,
Parece que agora funcionou. Obrigado. Não tenho certeza, mas me parece que
o atributo probability é a amplitude do intervalo de 95%, não é? Bom agora
estou aqui pensando ... com os dados origianis que levou quase dez minutos
reamostrando, talvez seja interessante um ajuste possivelmente com funções
do pacote paralelo como parSapply pra fazer mais rápido. Mas isso não deve
ser dificil.
> SMR2 <-
function(a){length(a$Hosp_Death[which(a$Hosp_Death==1)])/(length(a$SAPS3Pro2)
* mean(a$SAPS3Pro2))}
> smr.boot <- sapply(1:9999, function(x) SMR2(a[sample(1:nrow(a),
replace=TRUE),]))
> require(coda)
Carregando pacotes exigidos: coda
> hpdic <- HPDinterval(as.mcmc(smr.boot), 0.95)
>
> hist(smr.boot, breaks=pretty(smr.boot, 50))
> abline(v=hpdic)
> hpdic
lower upper
var1 0.9836823 1.016965
attr(,"Probability")
[1] 0.949995
>
Abraço forte,
Pedro Brasil
Em 23 de setembro de 2015 11:33, Elias Teixeira Krainski <
eliaskrainski@yahoo.com.br> escreveu:
> O argumento 'w' e' mudado internamente na funcao boot(), ou seja, nao e'
> usado como sendo igual a 1 internamente. Essa alteracao e' o truque usado
> para o processo de reamostragem considerando estratos.
>
> Creio que vc nao precisa usar o pacote boot() para fazer reamostragem. E
> pode fazer reamostragem simples, pois nao vi situacao de haver estratos nos
> seus dados. Assim, voce poderia fazer
>
> smr.fun <- function(x) sum(x$Hosp_Death)/sum(x$SAPS3Pro2)
> smr0 <- smr.fun(a)
> smr.boot <- sapply(1:9999, function(x) smr.fun(a[sample(1:nrow(a),
> replace=TRUE),]))
>
> ### HPD interval (Intervalo de mais alta densidade)
> require(coda)
> hpdic <- HPDinterval(as.mcmc(smr.boot), 0.95)
>
> hist(smr.boot, breaks=pretty(smr.boot, 50))
> abline(v=hpdic)
>
>
> Elias
>
> _______________________________________________
> R-br mailing list
> R-br@listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e fornea cdigo
> mnimo reproduzvel.
>
O pacote boot() usa paralelizar por default... essa solucao "manual" poderia substituir sapply() por parSapply() sim
Desculpem, eu ainda não consegui rodar por causa que ele não reconhece este objeto “a”. É necessário montar um arquivo e chamar com o read.table?
dput(a[sample(1:nrow(a),100),c("Hosp_Death","SAPS3Pro2")])
Error in dput(a[sample(1:nrow(a), 100), c("Hosp_Death", "SAPS3Pro2")]) : object 'a' not found --- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
Mauro, faz assim: ignore essa linha que inicia com dput, copie o que esta abaixo dessa linha, digite: a <- e cole o que voce copiou
Também estou com mesma dúvida que o Mauro, não consegui rodar e criar o objetivo a partir do dput() Enviado do Yahoo Mail no Android De:"Mauro Sznelwar" <sznelwar@uol.com.br> Data:21:34 Qua, 23 de Set de PM Assunto:[R-br] RES: Intervalos de confiança com boot.ci Desculpem, eu ainda não consegui rodar por causa que ele não reconhece este objeto “a”. É necessário montar um arquivo e chamar com o read.table?
dput(a[sample(1:nrow(a),100),c("Hosp_Death","SAPS3Pro2")])
Error in dput(a[sample(1:nrow(a), 100), c("Hosp_Death", "SAPS3Pro2")]) : object 'a' not found  Este email foi escaneado pelo Avast antivírus. www.avast.com _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne� c�igo m�imo reproduz�el.
participantes (5)
-
Alexandre Santos -
ASANTOS
-
Elias Teixeira Krainski -
Mauro Sznelwar -
Pedro Emmanuel Alvarenga Americano do Brasil