Sobre este documento

Este é um documento de finalidade didática que contém códigos úteis para pesquisas de dose eficaz. Para ter acesso a um ou mais documento (como codigos.R) enviar pedido para gguimaraes@unb.br ou helgabgs@unb.br.

library("Iso")
## Iso 0.0-17
library("boot")
library("ggplot2")
library("ez")
library(openxlsx)
source("C:\\Users\\gabri\\Downloads\\codigos.R")

Dados do estudo NCT02454868 para exemplo

bd=read.xlsx("C:\\Users\\gabri\\Downloads\\NCT02454868.xlsx")
bd=read.xlsx2(bd)

# g1 = pacientes em que houve sucesso, g2 = pacientes em que houve falha
g1=subset(bd,bd$qualidade.iot.adequada==TRUE)
g2=subset(bd,bd$qualidade.iot.adequada==FALSE)

# Amostra da variável qualidade.iot.adequada
bd$qualidade.iot.adequada
##  [1]  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
## [23]  TRUE FALSE FALSE  TRUE FALSE
# Amostra das doses
bd$dose
##  [1] 2.00 1.75 2.00 2.25 2.50 2.25 2.00 2.75 2.50 2.25 2.00 1.75 2.00 2.25
## [15] 2.50 2.75 3.00 3.25 3.00 3.25 3.00 2.75 3.00 2.75 3.00 3.25 3.00
# 
test.df=data.frame(responseSequence=bd$resp,doseSequence=bd$dose)

Código da figura da alocação sequencial

O código comentado abaixo foi usado para gerar a figura no formato eps, que é vetorial e comumente aceito por periódicos internacionais. O símbolo # foi usado para que ele não fosse executado nesta publicação e permitir que a figura seja mostrada no rpubs. Ainda não existe guia que determine qual cor deve ser usada para falhas ou sucesso, mas recomendamos preto e branco porque muitas revistas cobram taxa se figura colorida.

#setEPS()
#postscript("C:\\Users\\gabri\\Documents\\gabriel\\Figure2.eps")

plot(bd$dose~bd$id,data=bd,xlab="Patient number",ylab="Remifentanil mcg/kg")
lines(bd$dose~bd$id,data=bd)
points(g1$dose~g1$id,col="black",bg="white",pch=21) # bolas brancas para sucesso
points(g2$dose~g2$id,col="black",bg="black",pch=21) # bolas pretas para falhas

#dev.off()

Análise estatística com regressão isotônica e PAVA (código de exemplo)

No código abaixo, modifique a variável PercentilDE para o percentil adequado para o qual o estudo foi planejado. Por exemplo, para DE50 use PercentilDE=0.5. Para DE95 use PercentilDE=0.95. Importante: use percentis DE50 sempre que usar o método up-and-down de Dixon, para percentis diferentes existem vários métodos como 3+3, moeda enviesada (Biased Coin) e reavaliação contínua (continual reassessment method).

Não use regressões probit ou logísticas para esta finalidade porque há quebra de diversos pressupostos e prejuízo da análise, as melhores revistas não irão aceitar mais.

PercentilDE=0.5
testPava.df <- preparePava(test.df)
print(testPava.df)
##   naiveProbability pavaProbability nEvents nTrials nDoses
## 1        0.0000000       0.0000000       0       2   1.75
## 2        0.4000000       0.3333333       2       5   2.00
## 3        0.2500000       0.3333333       1       4   2.25
## 4        0.6666667       0.3846154       2       3   2.50
## 5        0.2500000       0.3846154       1       4   2.75
## 6        0.3333333       0.3846154       2       6   3.00
## 7        1.0000000       1.0000000       3       3   3.25
test.boot <- boot(data = test.df, statistic = bootIsotonicRegression, R = 999, sim = 'parametric', ran.gen = bootIsotonicResample, mle = list(baselinePava = testPava.df, firstDose = 2, PROBABILITY.GAMMA = 0.5), baselinePava = testPava.df, PROBABILITY.GAMMA = PercentilDE)
r2=bootBC.ci(test.boot$t0[3], test.boot$t[, 3])
ici=r2$`2.5% Bias Corrected Lower Bound`
DE=r2$`Original Statistic` #ED50
icm=r2$`97.5% Bias Corrected Upper Bound`
boot.ci(test.boot, type = c('norm', 'basic', 'perc'), conf = 0.95, index = 3)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = test.boot, conf = 0.95, type = c("norm", "basic", 
##     "perc"), index = 3)
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   ( 2.673,  3.759 )   ( 2.997,  4.021 )   ( 2.072,  3.097 )  
## Calculations and Intervals on Original Scale
print(paste("DE",PercentilDE*100,"(IC95%):",DE,"(",ici,"-",icm,")"))
## [1] "DE 50 (IC95%): 3.046875 ( 2.6304347826087 - 3.11363636363636 )"