6.6 Aplicações
1) Tabela 6.1 apresenta os resultados de um estudo de coorte, com vista a determinar, entre outras coisas, as vantagens de usar, como uma ferramenta de triagem, o teste ergométrico (EST) medidas. Em tais testes, um resultado passa/falha pode ser derivado durante o diagnóstico de doenças da artéria coronária (CAD) .
tabela 6.1. Screening test for physical stress
| CAD | ||||
|---|---|---|---|---|
| Healthy | Ill | Total | ||
| EST | Negative | 327 | 208 | 535 |
| Positive | 115 | 815 | 930 | |
| Total | 442 | 1 023 | 1 465 | |
Ele será assumido que não há verificação de viés:
a)
com Base nesta matriz de confusão, indicam os seguintes valores (com 95% de intervalo de confiança): sensibilidade e especificidade, o positivo e o valor preditivo negativo.
b)
Qual é o valor da área sob a curva para os dados atuais?
A criação da tabela de dados pode ser feito como indicado a seguir, a partir de uma tabela criada com o comando matrix() :
> tab < – como.tabela(matrix(c(815,115,208,327), nrow = 2, byrow = TRUE
dimnames = list(EST = c(“+”,”−”),
CAD = c(“+”,”-“))))
> tab
por Isso, deve ser observado que a tabela foi reorganizada, de modo a corresponder ao habitual notações, com “positivo” em eventos de primeira linha (normalmente, a exposição) e a primeira coluna (normalmente, a doença) da tabela.
seria então possível calcular todas as quantidades necessárias usando muito poucas operações aritméticas com R. No entanto, a embalagem epiR contém todos os instrumentos necessários para responder a questões epidemiológicas nos quadros 2 × 2. Assim, a epinefrina de comando.(testes) fornece prevalência, sensibilidade/especificidade, bem como os valores preditivos positivos e negativos:
> library(epiR)
> pav.testes (tab)
em relação à área sob a curva, um comando externo também pode ser empregado, roc.de.(tabela), no pacote epicalc:
> library(epicalc)
> roc.de.a tabela 6.2 (tab, graph = FALSE)
2) resume a proporção de enfartes do miocárdio observados em homens com idades compreendidas entre os 40 e os 59 anos e para os quais foram medidos o nível de pressão arterial e a taxa de colesterol, considerados sob a forma de classes ordenadas .
Tabela 6.2. Do miocárdio e da pressão arterial
| Colesterol (mg/100 ml) | |||||||
|---|---|---|---|---|---|---|---|
| BP | < 200 | 200 – 209 | 210 – 219 | 220 – 244 | 245 – 259 | 260 – 284 | > 284 |
| < 117 | 2/53 | 0/21 | 0/15 | 0/20 | 0/14 | 1/22 | 0/11 |
| 117 – 126 | 0/66 | 2/27 | 1/25 | 8/69 | 0/24 | 5/22 | 1/19 |
| 127 – 136 | 2/59 | 0/34 | 2/21 | 2/83 | 0/33 | 2/26 | 4/28 |
| 137 – 146 | 1/65 | 0/19 | 0/26 | 6/81 | 3/23 | 2/34 | 4/23 |
| 147 – 156 | 2/37 | 0/16 | 0/6 | 3/29 | 2/19 | 4/16 | 1/16 |
| 157 – 166 | 1/13 | 0/10 | 0/11 | 1/15 | 0/11 | 2/13 | 4/12 |
| 167 – 186 | 3/21 | 0/5 | 0/11 | 2/27 | 2/5 | 6/16 | 3/14 |
| > 186 | 1/5 | 0/1 | 3/6 | 1/10 | 1/7 | 1/7 | 1/7 |
The data are available in the file hdis.dat na forma de uma tabela com quatro colunas, indicando, respectivamente, a pressão arterial (oito categorias, com classificação de 1 a 8), a taxa de colesterol (sete categorias, com classificação de 1 a 7), o número de infarto do miocárdio e o número total de indivíduos. A associação entre a pressão arterial e a probabilidade de sofrer um infarto do miocárdio é o ponto de interesse:
a)
Calcular a proporção de infarto do miocárdio para cada nível de pressão arterial e representá-los em uma tabela e na forma gráfica;
b)
Expressar as proporções calculadas em (a), no logit formulário;
c)
com Base em um modelo de regressão logística, determinar se existe uma associação significativa no limiar α = 0.05 entre a pressão sangüínea, tratada como uma variável quantitativa, considerando a classe de centros e a probabilidade de ter um ataque cardíaco;
d)
Express no logit unidades de probabilidades de infarto do predito pelo modelo para cada um dos níveis de pressão arterial;
e)
Visualização no mesmo gráfico empírica proporções e a curva de regressão logística de acordo com os valores de pressão arterial (classe centros).importar os IDC dos dados.dat é alcançado como a seguir:
> bp< – read.quadro (“IDC.dat”, header = TRUE)
> str(bp)
Como pode ser verificado, não há nenhum rótulo associado com os níveis das variáveis de interesse (bpress para a pressão arterial, chol para a taxa de colesterol). Para gerar e associar etiquetas, os seguintes comandos podem ser inseridos:
> blab < – c(“<117”,”117-126”,”127-136”,”137-146”,
“147-156”,”157-166”,”167-186”,”>186”)
> clab < – c(“<200”,”200-209”,”210-219”,”220-244”,
“245-259”,”260-284”,”>284”)
> bp$bpress < – factor(bp$bpress, labels = blab)
> bp$chol < – factor(bp$chol, labels = clab)
O último comando converte as variáveis originais em variáveis qualitativas e associa às suas modalidades as etiquetas definidas por blab e clab. Para verificar que o banco de dados está agora na forma desejada, os seguintes comandos são úteis:
> str(bp)
> summary(bp)
A tabela de freqüência relativa expostas na declaração de agora pode ser reproduzido:
> round (xtabs (HDs/bpress total + chol, data = bp), 2)
Uma vez que vamos concentrar-nos apenas na relação entre pressão arterial e enfarte do miocárdio, é necessário agregar dados sobre os níveis de colesterol. Por outras palavras, é necessário acumular os valores para cada nível de pressão arterial, independentemente dos níveis de colesterol. A oportunidade também devem ser tomadas para mudar o nome os níveis da variável bpress usando os centros dos intervalos de classe:
> blab2 < – c(111.5,121.5,131.5,141.5,151.5,161.5,176.5,191.5)
> bp$bpress < rep(blab2, cada = 7)
> dfrm < agregada(bp
(lista de bpress = bp), soma)
é o último comando, de agregação (a), o que torna possível a agregação dos dados: – a soma de todos os valores (soma) de uma variável chol que não estão incluídos na lista de variáveis a ser retida para análise. Ao mesmo tempo, os resultados são armazenados em um novo banco de dados chamado dfrm. Note que neste caso não fazemos uso da notação da fórmula, mas da sintaxe alternativa para Agregado() (variável de resposta seguida por uma lista que consiste de uma ou mais variáveis explicativas). Uma visão geral dos dados agregados podem ser obtidos pela dando a cabeça() como entrada:
> head(dfrm, 5)
Quando uma proporção p está disponível, o seu valor em uma escala cujas unidades são logits é dada pela relação log(p/(1 − p)). Daí podem ser derivados os seguintes comandos para converter as proporções, calculadas como o IDIS / total em unidades logit:
> logit < – function(x) log(x/(1-x))
> dfrm$prop < – dfrm$idhs/dfrm$total
> dfrm$logit < – logit(dfrm$idhs/dfrm$total)
por Isso, deve ser observado que uma pequena função que tenha sido definido para a conversão de valores de x, que aqui é considerado um proporções, em seu equivalente log(x/(1 − x)). Da mesma forma, podia-se escrever:
> log((dfrm$idhs/dfrm$total)/(1-dfrm$idhs/dfrm$total))
Os valores calculados foram adicionados diretamente para o quadro de dados sob o nome de prop e logit.
O modelo de regressão logística é escrito da seguinte forma:
> glm(cbind(idh,o total de-idh) bpress,dados = dfrm,family = binomial)
A formulação a ser utilizada, cbind(idh, o total de-idh) bpress, leva em conta o fato de que estamos acessando dados agrupados e não respostas individuais. O comando glm() com a família de opções = binomial corresponde a uma regressão logística, que, sem entrar em detalhes, usa por padrão a função logit como um link canônico.
Os coeficientes de regressão podem ser apresentadas por meio da utilização do resumo():
> summary(glm(cbind(idh, o total de-idh) bpress, dados = dfrm,
family = binomial))
O resultado anterior inclui a informação essencial para responder à questão da significância estatística da associação entre pressão arterial e infarto do coração probabilidade: o coeficiente de regressão (na escala de log-odds) é igual a 0,024 e é significativo no limiar habitual de 5% (ver coluna Pr(>|z|)).
a probabilidade de ter um ataque cardíaco de acordo com os diferentes níveis de pressão arterial considerados é obtida da seguinte forma:
> glm.res < – glm(cbind(idh, o total de-idh) bpress, dados = dfrm,
family = binomial)
> predict(glm.res)
deve-se notar que os resultados intermédios gerados por R foram armazenados com o nome glm.res antes de usar o comando predict(). As previsões geradas por R são expressas em logit e é possível comparar os logits observados e previstos:
> cbind(dfrm, logit.predit = predict (glm.
praticamente todos os elementos estão disponíveis para representar as proporções de enfarte do miocárdio observadas e previstas graficamente de acordo com o nível de pressão arterial. Faltam as previsões do modelo expressas em proporções e não na escala logodds. Pode, portanto, ser desejável desenhar a curva de previsão, ou seja, a probabilidade de ter um ataque cardíaco de acordo com a pressão arterial, sem limitar esta última aos oito valores observados para o bpress variável. Aqui segue uma solução possível:
> dfrm$prop.predit < – predict (glm.res, type = “response”)
> f < – function(x)
1/(1 + exp(-(coeficiente(glm.res) + coef(glm.res)*x)))
> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,
xlab = “Blood pressure”,ylab = “Infarction likelihood”,
panel = function(x, y, …) {
panel.xyplot(x, y, col = “gray30”, pch = 19, …)
panel.curve(f, lty = 3, col = “gray70”)
panel.points(x, dfrm$prop.predit, col = “gray70”,…)
})
3) Uma investigação de caso-controlo centrou-se na relação entre o consumo de álcool e o consumo de tabaco e cancro esofágico em seres humanos (estudo”Ille e Villaine”). O grupo de casos consistiu em 200 pacientes que sofriam de câncer de esôfago, diagnosticados entre janeiro de 1972 e abril de 1974. Ao todo, 775 homens de controle foram selecionados das listas eleitorais. Quadro 6.3 mostra a distribuição da totalidade dos indivíduos de acordo com o seu consumo diário de álcool, considerando que um consumo superior a 80 g é considerado um fator de risco .
tabela 6.3. O consumo de álcool e câncer esofágico
| o consumo de Álcool (g/dia) | |||
|---|---|---|---|
| ≥ 80 | < 80 | Total | |
| Caso | 96 | 104 | 200 |
| Controle | 109 | 666 | 775 |
| Total | 205 | 770 | 975 |
a)
Qual é o valor de odds ratio e respectivo intervalo de confiança 95% (método de Woolf)? É uma boa estimativa do risco relativo?
b)
a proporção de consumidores em risco é a mesma entre os casos e entre os controlos (considerar α = 0,05)?
c)
construir o modelo de regressão logística tornando possível testar a associação entre o consumo de álcool e o Estatuto dos indivíduos. O coeficiente de regressão é significativo?
D)
encontrar o valor da razão de probabilidade observada, calculada em b) e o seu intervalo de confiança com base nos resultados da análise de regressão.uma vez que os dados brutos (individuais) não estão disponíveis, é necessário trabalhar diretamente com a tabela de frequência fornecida na declaração de problemas:
> alcohol < – matrix(c(666,104,109,96), nr = 2,
dimnames = list(c(“Control”,”Case”),
c(“<> = 80”)))
> alcohol
Regarding the odds ratio, the command oddsratio() from the package vcd will be used:
> library(vcd)
> oddsratio(álcool, log = FALSE)
A opção de log = FALSE garante que o resultado retornado corresponde corretamente ao odds ratio e não para o logaritmo da razão de probabilidades. Para obter uma assintótico intervalo de confiança a seguir, serão empregados:
> confint(oddsratio(álcool, log = FALSE))
da mesma forma, resumo(oddsratio(álcool)) pode ser usado para executar um teste de hipótese em que o log odds ratio (H0 : log(θ) = 0).
O suporte de comando.teste() pode ser usado para testar a hipótese de que a proporção de pessoas com um consumo diário é igual ou superior a 80 g e é idêntico nos casos, bem como nos controles, indicando os valores observados a partir da tabulação cruzada dada na declaração do problema:
> prop.teste(c(96,109), c(200,775), correto = FALSE)
Este teste é exatamente equivalente a um teste-Z para testar a diferença entre duas proporções estimadas a partir dos dados (se não a continuidade de correção está sendo empregado).
no que diz respeito ao modelo de regressão logística, a tabela de contingência tem de ser transformada numa tabela de dados em que as duas variáveis qualitativas (doença e exposição) estão representadas adequadamente, ou seja, com um quadro de dados. Isto pode ser feito usando o comando melt() do pacote reshape2 que torna possível transformar os dados em forma tabular em formato “longo”. Nesta fase, os níveis da variável de doença deverá ser alterada com 0 (controle) e 1 (casos), embora este não é realmente necessário e o nível 0 deve ser considerada como a categoria de referência (o que facilita a interpretação dos resultados):
> library(reshape2)
> álcool.df < – melt(álcool)
> nomes(álcool.df) < – c(“doença”, “exposição”, “n”)
> níveis(álcool.DF$disease) < – c(1,0)
> álcool.DF$disease < – relevel(alcohol.DF$disease, “0”)
O modelo de regressão logística é então escrito da seguinte forma:
> glm.res < – glm (exposição à doença, dados = álcool.DF,
família = binomial, weights = n)
> resumo(glm.res)
o resultado do interesse é a linha associada à exposição > = 80, uma vez que fornece informação sobre o valor do coeficiente de regressão associado à variável de exposição estimada por R, com o seu erro-padrão, bem como o valor da Estatística do ensaio. Aqui, o coeficiente de regressão é interpretado como o log da razão de probabilidades. Note – se que obteríamos exatamente os mesmos resultados trocando o papel das variáveis na formulação anterior, doença de exposição.
A razão de probabilidades calculada acima pode ser encontrada a partir do coeficiente de regressão associado ao factor de interesse (exposição), bem como do seu intervalo de confiança de 95%, da seguinte forma:
> exp(coef(glm.
> exp (confint (glm.res))
A segunda solução consiste em considerar o número de casos e o número total de indivíduos:
> alcohol2 < – dados.quadro(expos = c(“<> = 80”), caso = c(104,96),
total = c(104 + 666, 96 + 109))
> summary(glm(cbind(caso, o total de-caso) exposições, dados = alcohol2,
family = binomial))
4) estes dados vem de um estudo que procura estabelecer o prognóstico de validade de creatina quinase concentração no organismo na prevenção de ocorrência de infarto do miocárdio .
os dados estão disponíveis no ficheiro sck.dat: a primeira coluna corresponde à variável creatina cinase (ck), a segunda à presença da doença variável (pres) e a última à ausência variável de doença (abs).
A)
Qual é o número total de indivíduos?
b)
Compute the relative sick / healthy frequencies and represent their evolution according to creatine kinase values by using a scatterplot (points + segments connecting the points).
c)
com Base em um modelo de regressão logística, que visa prever a probabilidade de ficar doente, de calcular o valor de CK de que este modelo prevê que as pessoas sofrem da doença, considerando um limiar de 0,5 (se P (doente) ≥ 0.5, então doente = 1).
D)
graficamente representam as probabilidades de não serem previstas por este modelo, bem como as proporções empíricas de acordo com os valores ck.
E)
estabeleça a curva ROC correspondente e indique o valor da área sob a curva.
Desde que o nome das variáveis que não está no arquivo de dados, será necessário atribuir-lhes imediatamente depois que os dados são importados:
> sck < leitura.quadro (“sck.dat”, header = FALSE)
> nomes(sck) < – c(“ck”, “pres”, “abs”)
> summary(sck)
O número total de disciplinas corresponde à soma das frequências para as duas variáveis pres e abs, que é:
> sum(sck)
ou, equivalentemente: sum(sck$pres) + sum (sck$abs) (mas não sck$pres + sck$abs).
para calcular as frequências relativas destas duas variáveis, é necessário conhecer as quantidades totais por variável. Estes podem ser obtidos usando o comando aplicar e operacional por colunas:
> ni < – aplicar(sck, 2, sum)
com isso, basta dividir cada valor das variáveis pres e abs por os valores calculados anteriormente. Os valores obtidos serão armazenados em duas novas variáveis na mesma tabela de dados:
> sck$pres.prop < – sck$pres/ni
> sck$abs.prop < – sck$abs/ni
é possível verificar que os cálculos estão corretos: a soma dos valores para cada variável deve agora ser igual a 1:
> apply(sck, 2, sum)
em Seguida, as proporções obtidas simplesmente são representados em um único gráfico, considerando-se os valores da variável ck como o das abcissas:
> xyplot(pres.prop + abs.prop CK, data = sck, type = c (“b”,”g”),
auto.key = TRUE, ylab = “Frequency”)
the instruction type = c (“b”,” g”) means that we want the points to be displayed by lines (“b” = ” O ” + “l’) as well as a grid (“g”).
O modelo de regressão para dados agrupados é escrito como:
> glm.res < – glm(cbind(pres,abs) ck, dados = sck, family = binomial)
> summary(glm.res)
as previsões, expressas em forma de Probabilidades e não na escala de log-odds, são obtidas por meio do comando predict, especificando o tipo de opção = “resposta” da seguinte forma:
> glm.pred < – predict (glm.res, type = “response”)
> names(glm.pred) < – sck$ck
considerando que a probabilidade ≥ 0.5 designar pessoas que estão doentes, a seguinte distribuição é assim obtida:
> glm.pred
conclui-se que as pessoas serão consideradas doentes, de acordo com este modelo, para valores ck de 80 ou mais. Isso é facilmente verificado com um gráfico onde as probabilidades preditas com base nos valores da variável ck são representados:
> sck$doentes < – sck$pres/(sck$pres + sck$abs)
> xyplot(glm.pred ~ sck$ck, type = “l”,
ylab = “Probabilidade”, xlab = “ck”,
painel = função(…) {
painel.xyplot ( … )
painel.xyplot(sck$ck, sck$doentes, pch = 19, col = “cinza”)
})
Em um primeiro passo, é necessário “descompactar” os dados agrupados e para criar uma tabela com duas colunas: a primeira representa a variável de ck e o segundo representando a presença ou a ausência de doença. Usaremos as contagens por subgrupo previamente calculadas e disponíveis na variável ni:
> sck.expandir < – dados.quadro(ck = c(rep(sck$ck, sck$pres),
rep(sck$ck, sck$abs)),
doente = c(rep(1,ni), rep(0,ni)))
> tabela(sck.expanda$sick)
> com(sck.expanda, tapply (sick, ck, sum)
os dois últimos comandos são destinados a garantir que as mesmas contagens são encontradas e que a distribuição das pessoas doentes de acordo com o valor de ck é correta. Ele também pode ser verificado que os mesmos resultados são obtidos com relação a regressão logística e, ao mesmo tempo, os valores previstos na anterior tabela de dados podem ser adicionados (o mesmo procedimento anteriormente poderia ser seguido e as previsões podem ser replicadas, mas isso é mais simples desta forma):
> glm.res2 < – glm (sick ck, data = sck.expand, family = binomial)
> sck.expanda$prediction < – ifelse (predict (glm.res2,
Tipo = “resposta”) > = 0.5,
1, 0)
> com(sck.expanda, table (sick, prediction)
The last command makes it possible to display a confusion matrix in which the real diagnostics and those predicted by the regression model are crossed. As taxas de classificação corretas podem ser comparadas quando variamos o limiar de referência:
> classif.tab < – com (sck.expand, table(sick, prediction)
> sum(diag (classif.tab)) / sum (classif.tab)
Para exibir a curva ROC, vamos fazer uso do pacote ROCR:
> library(ROCR)
> pred < – predição(prever(glm.res2, type = “response”),
sck.expand$sick)
> perf < – performance(pred, “tpr”,”fpr”)
> plot(perf, ylab = “Sensitivity”, xlab = “1-Specificity”)
> grid()
> abline(0, 1, lty = 2)
The value of the area under the curve is obtained as follows:
> performance(pred, “auc”)@”y.values”