Capítulo 1 - Conceitos introdutórios

Exercício 1

Variável Resposta Variáveis Explicativas
(a) Infecção urinária Sexo e tratamento
(b) Câncer de esôfago Consumo de bebida alcoólica e histórico familiar
(c) Horas de alívio de dor de cabeça Dosagem do medicamento e idade
(d) Método de aprendizado preferido Período escolar

Exercício 2

Nominal Ordinal
Infecção urinária; Sexo; Tratamento Horas de alívio de dor de cabeça
Câncer de esôfago; Histórico familiar Dosagem do medicamento
Consumo de bebida alcoólica Idade
Método de aprendizado; Período escolar

Exercício 3

[a] estudo transversal
[b] entre pacientes do sexo masculino = 14/26 = 0,5385; e entre pacientes do sexo feminino = 12/13 = 0,9231

Exercício 4

[a] estudo de coorte
[b] entre fumantes = 90/800 = 0,1125; e entre não fumantes = 10/1200 = 0,0083

Exercício 5

[a] estudo caso-controle
[c] vantagens: custo e tempo reduzidos para obtenção da resposta quando comparado ao estudo de coorte.
desvantagens: propensos a vícios devido à possíveis manipulações dos grupos de comparação (casos e controles), bem como por estar condicionado à disponibilidade e qualidade dos registros existentes ou da memória dos pacientes.

Exercício 6

[a] ensaio clínico aleatorizado. Justificativa: a vitamina C foi administrada a um grupo de pacientes e ao outro não, caracterizando um estudo experimental e não observacional.

Exercício 7

[a] estudo caso-controle

Exercício 8

[a] estudo transversal

Capítulo 2 - Delineamentos amostrais e modelos associados

Exercício 1

[a] modelo multinomial

Exercício 2

Estudo descrito no exercício Modelo
4 Produto de binomiais
5 Produto de binomiais
6 Produto de binomiais
7 Produto de binomiais
8 Multinomial

Exercício 5

[b] Produto de binomiais

Exercício 6

[a] Produto de binomiais

Exercício 9

Sejam \(N_{11}\), \(N_{12}\), \(N_{21}\) e \(N_{22}\) varáveis aleatórias Poisson independentes com parâmetros \(\mu_{11}, \mu_{12}, \mu_{21}\) e \(\mu_{22}\), respectivamente. Então,

\(~~~ P(N_{11} = n_{11}, N_{12} = n_{12}, N_{21} = n_{21}, N_{22} = n_{22}) | N_{11}+N_{12}+N_{21}+N_{22}=n)\)=

= \(\displaystyle{\frac{P(N_{11}=n_{11},N_{12}=n_{12},N_{21}=n_{21},N_{22}= n -(n_{11}+n_{12}+n_{21}))}{P(N_{11}+N_{12}+N_{21}+N_{22}=n)}} = \frac{(1)}{(2)}\)

Como,

\((1) = \displaystyle{\frac{e^{-\mu_{11}}(\mu_{11})^{n_{11}}}{n_{11}!}\frac{e^{-\mu_{12}}(\mu_{12})^{n_{12}}}{n_{12}!} \frac{e^{-\mu_{21}}(\mu_{21})^{n_{21}}}{n_{21}!}\frac{e^{-\mu_{22}}(\mu_{22})^{n - (n_{11} + n_{12} + n_{21})}}{[n - (n_{11} + n_{12} + n_{21})]!}}\)

\((2) = \displaystyle{\frac{e^{-(\mu_{11}+ \mu_{12}+\mu_{21}+\mu_{22})}(\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22})^{n}}{n!}}\)

segue que,

\(\displaystyle{\frac{(1)}{(2)} = \frac{n!}{n_{11}!~n_{12}!~n_{21}!~[n - (n_{11}+n_{12}+n_{21})]!}~ \frac{(\mu_{11})^{n_{11}}(\mu_{12})^{n_{12}}(\mu_{21})^{n_{21}}(\mu_{22})^{n - (n_{11}+n_{12}n_{21})}} {(\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22})^n}}\)

\(= \frac{n!}{n_{11}!n_{12}!n_{21}![n - (n_{11}+n_{12}+n_{21})]!}~(\frac{\mu_{11}}{\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22}})^{n_{11}}~(\frac{\mu_{12}}{(\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22})})^{n_{12}}~(\frac{\mu_{21}}{(\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22})})^{n_{21}} ~(\frac{\mu_{22}}{(\mu_{11}+\mu_{12}+\mu_{21}+\mu_{22})})^{n -(n_{11}+n_{12}+n_{21})}\)

\(= \displaystyle{\frac{n!}{n_{11}!~n_{12}!~n_{21}!~[n - (n_{11}+n_{12}+n_{21})]!} (p_{11})^{n_{11}} (p_{12})^{n_{12}} (p_{21})^{n_{21}} (p_{22})^{n - (n_{11}+n_{12}+n_{21})}}\)

em que \(p_{ij} = \frac{\mu_{ij}}{\sum_{i,j = 1}^2 \mu_{ij}}\) para i, j = 1,2.

Capítulo 3 - Tabelas de contingência 2 x 2

Exercício 1

1.1 Dados do exercício 3 - estudo transversal    
[a] Proporções amostrais e frequências esperadas

dados<-matrix(c(14,12,12,1), nc=2)
n <- sum(dados)
dados/n
          [,1]       [,2]
[1,] 0.3589744 0.30769231
[2,] 0.3076923 0.02564103

Qp<-chisq.test(dados)
Qp$expected 
          [,1]     [,2]
[1,] 17.333333 8.666667
[2,]  8.666667 4.333333

[b] Hipótese nula de interesse => H0: hipótese de independência
 
fisher.test(dados, alternative="two.sided")   # bilateral
p-value = 0.02857
 
fisher.test(dados, alternative="less")        # unilateral
p-value = 0.01674
 
[c] Condicional aos totais marginais-linha
        OR         li       ls
  0.09722222 0.01098085 0.860786

Exercício 2

[a] Estudo de coorte
[c]    RR        li       ls
    0.9440171 0.6826663 1.305423

Exercício 3

[b] Qui-quadrado de Pearson e RR  

dados<-matrix(c(83,34,19,16), nc=2)
Qp<-chisq.test(dados, correct=F) 
Qp  

X-squared = 3.3852, df = 1, p-value = 0.06578
IC(RR)_95%        RR        li       ls
            1.306793 0.9438295 1.809338
IC(RR)_90%        RR        li       ls
            1.306793 0.9945187 1.717119

Exercício 4

Para \(g_1(\widehat \theta) = \ln(\widehat p_{(1)1})\) e \(g_2(\widehat \theta) = \ln(\widehat p_{(2)1})\), segue pelo método delta que

\(Var[g_1(\widehat \theta)] = Var[\ln(\widehat p_{(1)1}] = \displaystyle{\frac{p_{(1)1}(1-p_{(1)1})}{n_{1+}}\bigg(\frac{1}{p_{(1)1}}\bigg)^2 = \frac{(1-p_{(1)1})}{(n_{1+})p_{(1)1}}}\)

\(Var[g_2(\widehat \theta)] = Var[\ln(\widehat p_{(2)1}] = \displaystyle{\frac{p_{(2)1}(1-p_{(2)1})}{n_{2+}}\bigg(\frac{1}{p_{(2)1}}\bigg)^2 = \frac{(1-p_{(2)1})}{(n_{2+})p_{(2)1}}}\)

Logo,

\(Var(\widehat f) = Var[\ln(\widehat p_{(1)1}] + Var[\ln(\widehat p_{(2)1}] = \displaystyle{ \frac{(1-p_{(1)1})}{(n_{1+})p_{(1)1}} + \frac{(1-p_{(2)1})}{(n_{2+})p_{(2)1}}}\)

Exercício 6

Terapia Óbito (Sim) Óbito (não) Total
Medicamentosa 42 161 203
Cirurgia 22 181 203
Angioplastia 29 176 205 
dados<-matrix(c(42,22,29,161,181,176), nc=2)
Qp<-chisq.test(dados, correct=F) 
Qp

X-squared = 7.911, df = 2, p-value = 0.01915
                                               OR      li       ls
Medicamentosa x Cirurgia        IC(OR)_95%   2.14624  1.22863   3.74917  
Medicamentosa x Angioplastia    IC(OR)_90%   1.58320  1.02397   2.44786  
Cirurgia x Angioplastia         IC(OR)_95%   0.73766  0.40820   1.33302

Exercício 7

X-squared = 1.5611, df = 1, p-value = 0.2115
 
           OR        li      ls
[1,] 1.133493 0.9311736 1.37977

Exercício 8

X-squared = 108.67, df = 2, p-value < 2.2e-16
           
                                     OR        li        ls
Usual x restrição sal              0.13063   0.08631   0.19772 
Usual x restrição gordura          0.32189   0.22479   0.46094  
Restrição sal x restrição gordura  2.46408   1.61405   3.76179   

                                     RR        li        ls
Usual x restrição sal              0.49175   0.42190   0.57316
Usual x restrição gordura          0.60356   0.51147   0.71224  
Restrição sal x restrição gordura  1.22737   1.11458   1.35158

Capítulo 4 - Tabelas de contingência s x r

Exercício 1

[a] \(H_0: \bar F_1 = \bar F_2\) versus \(H_A: \bar F_1 \neq \bar F_2\)

     fb1      fb2       Qs  gl           p
 3.086207 2.086957 12.21358  1    0.0004744286

[b] \(H_0: \bar F_1 = \bar F_2\) versus \(H_A: \bar F_1 \neq \bar F_2\)

     fb1      fb2       Qs  gl           p
 1.543103 1.043478 12.21358  1    0.0004744286

Exercício 2

\(H_0: p_{(1)1} = p_{(2)1} = p_{(3)1} = p_{(4)1} = p_{(5)1}\) (hipótese de homogeneidade)

\(H_A: p_{(i)1} \neq p_{(i')1}\) para pelo menos um par \((i, i')\) com \(i \neq i'\) (\(i, i'= 1,2,3,4,5\))

dados<-matrix(c(26,26,23,18,9,6,7,9,14,25), nc=2) 
Qp<-chisq.test(dados,correct=F) 
Qp 
Pearson's Chi-squared test: X-squared = 29.1235, df = 4, p-value = 7.378e-06

Para proceder as comparações “duas a duas” pode-se obter e comparar

\((a)\) IC_95% de \(p_{(i)1}\) ou de \(p_{(i)2}\), para \(i = 1,2,3,4,5\)
\((b)\) IC_95% de RR i|i’, com \(i \neq i'\) (\(i, i'= 1,2,3,4,5\))
\((c)\) IC_95% de OR i|i’, com \(i \neq i'\) (\(i, i'= 1,2,3,4,5\))

Utilizando (a) segue que:

\(\hat p_{(1)2}\) = 0,1875 com IC = (0,05; 0,32)
\(\hat p_{(2)2}\) = 0,2121 com IC = (0,07; 0,35)
\(\hat p_{(3)2}\) = 0,28125 com IC = (0,12; 0,44)
\(\hat p_{(4)2}\) = 0,4375 com IC = (0,26; 0,61)
\(\hat p_{(5)2}\) = 0,7353 com IC = (0,58; 0,88)

Ainda, assumindo o vetor de escores a = (0,1) e, se fizer sentido assumir o vetor de escores c = (0,1,2,3,4), pode-se testar

\(H_0\): ausência de tendência linear (rac = 0) versus \(H_A\): presença de tendência linear (rac\(\neq\) 0)

x<-c(rep(0,32),rep(1,33),rep(2,32),rep(3,32),rep(4,34))  
y<-c(rep(0,26),rep(1,6),rep(0,26),rep(1,7),rep(0,23),rep(1,9), rep(0,18),rep(1,14),rep(0,9),rep(1,25)) 
rac<-cor(y,x)  
n<-length(x) 
QCS<-(n-1)*rac^2  
p<-1-pchisq(QCS,1)
cbind(rac, QCS, p)

       rac      QCS            p
  0.391568 24.83873 6.233236e-07

Exercício 3

\(H_0: \bar F_1 = \bar F_2\) versus \(H_A: \bar F_1 \neq \bar F_2\)

dados<-matrix(c(4,6,4,8,4,2,6,4),nc=4)
escore<-c(0,1,2,3)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3]),sum(dados[,4])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs = ((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1
p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p)
 
 fb1      fb2       Qs gl         p
 1.666667 1.2 1.549573  1 0.2131985

Exercício 4

\(H_0: r_{ac} = 0\) versus \(H_A: r_{ac} \neq 0\)

x<-c(rep(6,21),rep(8.5,49),rep(10.5,50),rep(12.5,59),rep(14.5,44))  
y<-c(rep(1,7), rep(2,4), rep(3,3), rep(4,7), 
     rep(1,10),rep(2,15),rep(3,11),rep(4,13),
     rep(1,23),rep(2,9), rep(3,11),rep(4,7),
     rep(1,28),rep(2,9), rep(3,12),rep(4,10),
     rep(1,32),rep(2,5), rep(3,4), rep(4,3)) 
rac<-cor(y,x)  
n<-length(x) 
QCS<-(n-1)*rac^2  
p<-1-pchisq(QCS,1)
cbind(rac, QCS, p)

         rac      QCS            p
  -0.2801432 17.42261 2.992459e-05

Exercício 5

\(H_0: \bar F_1 = \bar F_2\) versus \(H_A: \bar F_1 \neq \bar F_2\)

    fb1      fb2            Qs   gl  p
    1.747097 1.258457 72.24612   1   0

Exercício 6

A versus B        fbA   fbB          Qs   gl        valor p
                  6.25  3.5    47.38783    1    5.824341e-12
A versus C        fbA   fbC          Qs   gl        valor p
                  6.25  4.86   21.07151    1    4.424585e-06
A versus D        fbA   fbD          Qs   gl        valor p
                  6.25  7.42   17.04376    1    3.652824e-05
B versus C        fbB   fbC          Qs   gl        valor p
                  3.5   4.86   18.52317    1    1.678513e-05
B versus D        fbB   fbD          Qs   gl        valor p
                  3.5   7.42   67.80638    1    2.220446e-16
C versus D        fbC   fbD          Qs   gl        valor p
                  4.86  7.42   51.39975    1    7.535084e-13

Capítulo 5 - Análise estratificada

Exercício 1

 dados<-matrix(c(45,80,84,106,64,104,124,117,71,116,82,87),nc=3)
 dados
 escore<-c(0,1,2) 
 fb11<-(sum(dados[1,]*escore))/sum(dados[1,]) 
 fb12<-(sum(dados[2,]*escore))/sum(dados[2,]) 
 fb21<-(sum(dados[3,]*escore))/sum(dados[3,]) 
 fb22<-(sum(dados[4,]*escore))/sum(dados[4,]) 
 fm1<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(fb11,fb21)) 
 esp1<-(c(sum(dados[1:2,1]),sum(dados[1:2,2]),sum(dados[1:2,3])))/sum(dados[1:2,]) 
 mu1<-sum(escore*esp1) 
 esp2<-(c(sum(dados[3:4,1]),sum(dados[3:4,2]),sum(dados[3:4,3])))/sum(dados[3:4,]) 
 mu2<-sum(escore*esp2) 
 mu<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(mu1,mu2)) 
 v1<- sum(((escore-mu1)^2)*esp1) 
 v2<- sum(((escore-mu2)^2)*esp2) 
 vfma<-(sum(dados[1,])*sum(dados[2,])*v1)/(sum(dados[1:2,])-1) 
 vfmb<-(sum(dados[3,])*sum(dados[4,])*v2)/(sum(dados[3:4,])-1) 
 vfm<- sum(c(vfma,vfmb)) 
 QSMH<-((fm1-mu)^2)/vfm 
 p<-1-pchisq(QSMH,1) 
 round(c(QSMH,p),digits=5)
 
 QSMH       p
 0.73789 0.39034

Exercício 2

 tab<-array(c(21,52,538,636,47,123,1578,1571),dim=c(2,2,2))
 mantelhaen.test(tab, correct=F)
 
 Mantel-Haenszel chi-squared test without continuity correction
 Mantel-Haenszel X-squared = 40.117, df = 1, p-value = 2.392e-10
 alternative hypothesis: true common odds ratio is not equal to 1
 95 percent confidence interval:
  0.3059031 0.5423354
 
 sample estimates:
 common odds ratio
         0.4073107 
 source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/breslowday.test.R.txt")
 breslowday.test(tab)

 Breslow and Day test (with Tarone correction):
 Breslow-Day X-squared         = 0.5116056
 Breslow-Day-Tarone X-squared  = 0.5114644
 Test for test of a common OR: p-value =  0.4745056
Fem|Masc: ORsem = 0.288  e  OR1 = 0.477  OR2 = 0.380   ORcomum = 0.407
Masc|Fem: ORsem = 3.466  e  OR1 = 2.095  OR2 = 2.628   ORcomum = 2.455

Exercício 3

a) Pacientes sob diagnóstico I

           rac     QCS              p
       0.391568   24.83873 6.233236e-07
(b) Pacientes sob diagnóstico II
 
           rac     QCS              p
      0.5718972   52.33063 4.689582e-13
(c) Controlando por diagnóstico  
  
    QCSMH = 67.14238 p < 0.0001

Os resultados evidenciam associação entre efeitos adversos e dosagens do medicamento. Há um aumento na proporção de pacientes com efeitos adversos (rac > 0) com o aumento da dosagem do medicamento. Ainda, rac(I+II) = 0,456, rac(I) = 0,3915 e rac(II) = 0,5718, indicando que a proporção de pacientes com efeitos adversos depende do diagnóstico.

x<-c(rep(0,64),rep(1,65),rep(2,65),rep(3,64),rep(4,66))
y<-c(rep(0,52),rep(1,12),rep(0,38),rep(1,27),rep(0,36),rep(1,29),rep(0,19),rep(1,45),rep(0,10),rep(1,56)) 
rac<-cor(y,x) 
n<-length(x) 
QCS<-(n-1)*rac^2 
p<-1-pchisq(QCS,1) 
cbind(rac,QCS,p) 

       rac      QCS            p
[1,] 0.4559289 67.14238 2.220446e-16

Exercício 4

(a) Associação entre fumo voluntário e câncer de pulmão, controlando por fumo passivo
 
tab<-array(c(120,111,80,155,161,117,130,124),dim=c(2,2,2))
mantelhaen.test(tab, correct=F)
 
Mantel-Haenszel chi-squared test without continuity correction
Mantel-Haenszel X-squared = 14.4223, df = 1, p-value = 0.0001461
 
Obs: Há evidências de associação entre fumo voluntário e câncer de pulmão, controlando por fumo passivo.
 
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/breslowday.test.R.txt")
breslowday.test(tab)

Breslow-Day X-squared                 = 3.273719
Test for test of a common OR: p-value = 0.07040859
ORsem = 1.63    OR1 = 2.09   OR2 = 1.31

Obs: a razão de chances comum não é apropriada para esse estudo, pois as razões de chances individuais sugerem chance de câncer de pulmão superior entre os fumantes voluntários que também são fumantes passivos, do que entre os fumantes voluntários que não são fumantes passivos.

Capítulo 6 - Tabelas com dados relacionados

Exercício 1

 dados<-matrix(c(40,8,10,70), nc=2)
 mcnemar.test(dados, correct=F)

 McNemar's chi-squared = 0.22222, df = 1, p-value = 0.6374

Exercício 2

 dados<-matrix(c(5,50,30,4), nc=2)
 mcnemar.test(dados, correct=F)

 McNemar's chi-squared = 5, df = 1, p-value = 0.02535

Exercício 3

 dados<-matrix(c(10,5,98,7), nc=2)
 mcnemar.test(dados, correct=F)

 McNemar's chi-squared = 83.971, df = 1, p-value < 2.2e-16

Exercício 5

 (a) Concordância quanto à gravidade: residente versus especializando
 
 require(vcd)              # instalar pacote vcd
 x<-c(89,5,1,1,3,0,0,1,0) 
 x<-matrix(x,3,3) 
 Kappa(x)
             value    ASE     z Pr(>|z|)
 Unweighted 0.4338 0.1506 2.880 0.003975
 Weighted   0.4368 0.1456 3.001 0.002691
(b) Concordância quanto à resolução: residente versus especializando
 
 require(vcd)
 x<-c(10,1,36,0,1,1,0,1,50) 
 x<-matrix(x,3,3) 
 Kappa(x)
             value     ASE     z  Pr(>|z|)
 Unweighted 0.2333 0.06439 3.624 0.0002900
 Weighted   0.2257 0.06349 3.555 0.0003782              
(c) Concordância quanto à gravidade da doença: mães versus residente
 
 require(vcd)
 x<-c(1,1,8,0,1,2,5,6,73)  
 x<-matrix(x,3,3)
 Kappa(x)
              value    ASE      z Pr(>|z|)
 Unweighted 0.09307 0.1003 0.9280   0.3534
 Weighted   0.07920 0.1053 0.7521   0.4520

Exercício 6

(a)  Concordância entre os examinadores -  Adesivo A

|              |  Examinador 2 |   |              | Examinador 3  |   |              |  Examinador 3 |
| Examinador 1 | 0  1  2  3  4 |   | Examinador 1 | 0  1  2  3  4 |   | Examinador 2 | 0  1  2  3  4 |
|------------------------------|   |------------------------------|   |------------------------------|
|     0        | 3  1  0  0  0 |   |     0        | 3  1  0  0  0 |   |     0        | 3  0  0  0  0 |
|     1        | 0  5  0  0  0 |   |     1        | 0  4  1  0  0 |   |     1        | 0  5  1  0  0 |
|     2        | 0  0  0  1  0 |   |     2        | 0  0  1  0  0 |   |     2        | 0  0  0  0  0 |
|     3        | 0  0  0  0  0 |   |     3        | 0  0  0  0  0 |   |     3        | 0  0  1  0  0 |
|     4        | 0  0  0  0  4 |   |     4        | 0  1  0  0  3 |   |     4        | 0  1  0  0  3 |

 Examinador 1 versus Examinador 2 - Adesivo A  =>  kw = 0,9152
 Examinador 1 versus Examinador 3 - Adesivo A  =>  kw = 0,7756
 Examinador 2 versus Examinador 3 - Adesivo A  =>  kw = 0,7756  
 
 require(vcd)
 x<-c(3,0,0,0,0,1,5,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,4)
 x<-matrix(x,5,5)
 summary(Kappa(x))

 x<-c(3,0,0,0,0,1,4,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,3)
 x<-matrix(x,5,5) 
 summary(Kappa(x)) 

 x<-c(3,0,0,0,0,0,5,0,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,3) 
 x<-matrix(x,5,5) 
 summary(Kappa(x))
 Concordância entre os examinadores -  Adesivo B

|              |  Examinador 2 |   |              | Examinador 3  |   |              |  Examinador 3 |
| Examinador 1 | 0  1  2  3  4 |   | Examinador 1 | 0  1  2  3  4 |   | Examinador 2 | 0  1  2  3  4 |
|------------------------------|   |------------------------------|   |------------------------------|
|     0        | 2  3  0  0  0 |   |     0        | 3  2  0  0  0 |   |     0        | 2  0  0  0  0 |
|     1        | 0  1  1  0  2 |   |     1        | 0  2  2  0  0 |   |     1        | 1  2  1  1  0 |
|     2        | 0  0  0  0  0 |   |     2        | 0  0  0  0  0 |   |     2        | 0  0  2  0  1 |
|     3        | 0  1  0  0  2 |   |     3        | 0  0  0  2  0 |   |     3        | 0  0  0  0  0 |
|     4        | 0  0  2  0  0 |   |     4        | 0  0  1  0  1 |   |     4        | 0  2  0  1  1 |

 Examinador 1 versus Examinador 2 - Adesivo B  =>  kw = 0,236
 Examinador 1 versus Examinador 3 - Adesivo B  =>  kw = 0,688
 Examinador 2 versus Examinador 3 - Adesivo B  =>  kw = 0,401
(b) Grau de infiltração entre os adesivos A e B

|              |      Grau     |   |              |     Grau      |   |              |      Grau     |
| Examinador 1 | 0  1  2  3  4 |   | Examinador 2 | 0  1  2  3  4 |   | Examinador 3 | 0  1  2  3  4 |
|------------------------------|   |------------------------------|   |------------------------------|
| Adesivo A    | 4  5  1  0  4 |   | Adesivo A    | 3  6  0  1  4 |   | Adesivo A    | 3  6  2  0  3 |
| Adesivo B    | 5  4  0  3  2 |   | Adesivo B    | 2  5  3  0  4 |   | Adesivo B    | 3  4  3  2  2 |

\(H_0: \bar F_1 = \bar F_2\) versus \(H_A: \bar F_1 \neq \bar F_2\)

Examinador 1 => \(\bar f_A = 1,64; \bar f_B =1,50; Q_S\) = 0,0577 (g.l. = 1, valor p = 0,8102)

Examinador 2 => \(\bar f_A = 1,79; \bar f_B =1,93; Q_S\) = 0,0608 (g.l. = 1, valor p = 0,8052)

Examinador 3 => \(\bar f_A = 1,57; \bar f_B =1,71; Q_S\) = 0,0735 (g.l. = 1, valor p = 0,7862)

Exercício 8

dados<-matrix(c(7,19,7,128),nc=2)
mcnemar.test(dados,correct=F)
 
    McNemar's Chi-squared test
McNemar's chi-squared = 5.5385, df = 1, p-value = 0.0186

Exercício 9

(a) Código para obtenção do gráfico de concordância 

require(vcd)
x<-c(22,5,0,0,0, 2,7,2,1,0, 2,14,36,14,3, 0,0,0,7,0, 0,0,0,0,3)
x<-matrix(x,5,5)
rownames(x) <- c("1","2","3","4","5")
colnames(x) <- c("1","2","3","4","5") 
x
agreementplot(x, xlab ="Categorias: patologista 1", ylab ="Categorias: patologista 2", line_col="gray")
   
(b) Intensidade de concordância

require(vcd)
x<-c(22,5,0,0,0, 2,7,2,1,0, 2,14,36,14,3, 0,0,0,7,0, 0,0,0,0,3)
x<-matrix(x,5,5)
Kappa(x)
summary(Kappa(x))
confint(Kappa(x), level=0.95)

            value     ASE      z  Pr(>|z|)
Weighted   0.6492 0.04867 13.339 1.369e-40

 Kappa           lwr       upr 
Weighted   0.5538055 0.7445806

Exercício 10

(b) Concordância entre os dois patologistas

require(vcd)
x<-c(291,186,2,3,1, 74,256,4,10,7, 1,7,0,1,1, 1,7,2,14,8, 1,3,0,2,3)
x<-matrix(x,5,5)
Kappa(x)
summary(Kappa(x))
confint(Kappa(x), level=0.95)

            value     ASE     z  Pr(>|z|)
 Weighted   0.4359 0.02834 15.38 2.074e-53
 
  Kappa           lwr       upr
 Weighted   0.3804014 0.4914762

Capítulo 7 - Regressão Binomial

Exercício 1

(a) Resultados do modelo de regressão logística

chd<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex1cap7.txt", h=T)
attach(chd)
mod1<-glm(as.matrix(chd[,c(1,2)]) ~ cat+idade+ecg+cat*idade+cat*ecg+idade*ecg, family=binomial(link="logit"),           data=chd)
anova(mod1, test="Chisq")
 
          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)   
NULL                          7    21.3320             
cat        1  14.1312         6     7.2008 0.0001705 ***
idade      1   4.7240         5     2.4768 0.0297448 * 
ecg        1   1.5224         4     0.9544 0.2172564   
cat:idade  1   0.0321         3     0.9223 0.8577925   
cat:ecg    1   0.5035         2     0.4188 0.4779635   
idade:ecg  1   0.4157         1     0.0031 0.5190999
 
Conclusão sobre o efeito das interações duplas => não significativos

# Modelo sem as interações (modelo 2) mostrou efeito não significativo de ecg na presença de cat e idade

mod2<-glm(as.matrix(chd[,c(1,2)]) ~ cat+idade+ecg, family=binomial(link="logit"), data=chd)
anova(mod2, test="Chisq")
summary(mod2)

# Modelo selecionado (modelo 3) 

mod3<-glm(as.matrix(chd[,c(1,2)]) ~ cat+idade, family=binomial(link="logit"), data=chd)
anova(mod3, test="Chisq")
summary(mod3)
       
dev<-residuals(mod3,type='deviance')
dev
QL<-sum(dev^2)
p1<-1-pchisq(QL,5)
cbind(QL,p1)

rpears<-residuals(mod3,type='pearson')
rpears
QP<-sum(rpears^2)
p2<-1-pchisq(QP,5)
cbind(QP,p2)

par(mfrow=c(1,2))
plot(rpears, ylab="resíduos Pearson", pch=16, ylim=c(-2,2))
abline(h=0,lty=3)   
plot(dev, ylab="resíduos deviance", pch=16, ylim=c(-2,2))
abline(h=0,lty=3)   
 
par(mfrow=c(1,1))
ntot<-c(274,122,59,32,8,39,17,58)
fit.model<-mod3
source("http://www.ime.usp.br/~giapaula/envelr_bino")
# Conclusões com base no modelo 3

 Coefficients: mod3
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.5397     0.2005 -12.664  < 2e-16 ***
cat           0.7736     0.2955   2.618  0.00884 ** 
idade         0.6174     0.2840   2.173  0.02975 *  

> exp(0.7736)        # fixando idade e variando o nível de catecholamine
[1] 2.167555
 
Pacientes com cat = 1 (nível de catecholamine alto) apresentaram chance de doença coronária igual a 2,16 vezes a dos pacien- tes com cat = 0 (nível de catecholamine baixo).  
 
> exp(0.6174)         # fixando nível de catecholamine e variando a idade
[1] 1.854101     
 
Pacientes com idade = 1 (>= 55 anos) apresentaram chance de doença coronária igual a 1,85 vezes a dos pacientes com idade = 0 (< 55 anos).  
 
> exp(0.7736 + 0.6174)  # variando a idade e o nível de catecholamine
[1] 4.018867
 
Pacientes com cat = 1 (nível de catecholamine alto) e idade = 1 (>= 55 anos) apresentaram chance de doença coronária igual a 4 vezes a dos pacientes com cat = 0 (nível de catecholamine baixo) e idade = 0 (< 55 anos).

Exercício 3

(a) script para obtenção das probabilidades ordenadas

 x1<-rep(0:5,13)
 x2<-rep(3:15,rep(6,13))
 px<-(exp(2.211 + 2.607*x1 - 0.52*x2))/(1+ exp(2.211 + 2.607*x1 - 0.52*x2))
 est<-as.data.frame(cbind(x1,x2,px))
 i<-order(est[,3], decreasing=T)
 px_ord<-est[i,]
 px_ord

Exercício 4

# Selecionando o modelo (considerado log10(dose))
 
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex4cap7.txt", h=T)
dados$ldose<-log10(dados$dose)
attach(dados)
 
ajuste1<-glm(as.matrix(dados[,c(3,4)])~droga+ldose+droga*ldose,family=binomial(link="logit"),data=dados)
anova(ajuste1,test="Chisq")
ajuste1<-glm(as.matrix(dados[,c(3,4)])~ldose,family=binomial(link="logit"),data=dados)
anova(ajuste1,test="Chisq")
 
ajuste2<-glm(as.matrix(dados[,c(3,4)])~droga+ldose+droga*ldose,family=binomial(link="cloglog"),data=dados)
anova(ajuste2,test="Chisq")
ajuste2<-glm(as.matrix(dados[,c(3,4)])~ldose,family=binomial(link="cloglog"),data=dados)
anova(ajuste2,test="Chisq")
summary(ajuste2)
 
ajuste3<-glm(as.matrix(dados[,c(3,4)])~droga+ldose+droga*ldose,family=binomial(link="probit"),data=dados)
anova(ajuste3,test="Chisq")
ajuste3<-glm(as.matrix(dados[,c(3,4)])~ldose,family=binomial(link="probit"),data=dados)
anova(ajuste3,test="Chisq")
summary(ajuste3)
 
ajuste4<-glm(as.matrix(dados[,c(3,4)])~droga+ldose+droga*ldose,family=binomial(link="cauchit"),data=dados)
anova(ajuste4,test="Chisq")
ajuste4<-glm(as.matrix(dados[,c(3,4)])~ldose,family=binomial(link="cauchit"),data=dados)
anova(ajuste4,test="Chisq")
summary(ajuste4)
 
# Obs: efeito da interação droga*ldose e efeito da droga não significativos
 
cbind(ajuste1$aic, ajuste2$aic, ajuste3$aic, ajuste4$aic)

       [,1]    [,2]     [,3]     [,4]
[1,] 29.75466 32.3345 29.68889 31.56915
 
# Obs: modelo binomial com link probito e o com link logito apresentaram os menores AIC  
 
par(mfrow=c(1,4))
ntot<-c(20,20,20,20,20,20,20,20)
fit.model<-ajuste1
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste2
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste3
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste4
source("http://www.ime.usp.br/~giapaula/envelr_bino")

# Modelos plausíveis: probito e logito  
Nota: se for considerando o modelo probito (ajuste 3) com x = log10(dose), tem-se:
    
summary(ajuste3)    
Deviance Residuals:
     Min        1Q    Median        3Q       Max 
-0.97037  -0.49554   0.08746   0.24303   0.99555 
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -2.1924     0.3183  -6.889 5.63e-12 ***
ldose         2.9665     0.3975   7.462 8.52e-14 ***
 
    Null deviance: 75.0849  on 7  degrees of freedom
Residual deviance:  3.0751  on 6  degrees of freedom
AIC: 29.689
 
ajuste3$y
   1    2    3    4    5    6    7    8
0.10 0.45 0.70 0.95 0.05 0.30 0.70 0.85
ajuste3$fitted.values         
1          2          3          4          5          6          7          8
0.09690486 0.34223190 0.68674050 0.91615184 0.09690486 0.34223190 0.68674050 0.91615184
ldose
[1] 0.30103 0.60206 0.90309 1.20412 0.30103 0.60206 0.90309 1.20412
dose
[1]  2  4  8 16  2  4  8 16
 
     Por exemplo, os camundongos que receberam dosagem 16 (log10(dose) = 1.20412) apresentaram chance de morte igual a aproximadamente 5 vezes a daqueles que receberam dosagem 8 (log10(dose) = 0.90309).

     OR 16|8 = [0.91615184/(1 - 0.91615184)]/[0.68674050/(1 - 0.68674050)] = 4.984086

----------------------------------------------------------------------- 

NOTA: se for considerado o modelo logito (ajuste 1) com x = log10(dose), 
tem-se OR 16|8 = exp[5.0357*(1.20412-0.90309)] = 4.553503.

summary(ajuste1)
Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.92315  -0.47382   0.00253   0.16679   1.04263  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.7070     0.5925  -6.256 3.94e-10 ***
ldose         5.0357     0.7511   6.704 2.03e-11 ***

    Null deviance: 75.0849  on 7  degrees of freedom
Residual deviance:  3.1409  on 6  degrees of freedom
AIC: 29.755

Exercício 6

(a) Ajustando e selecionando um modelo
 
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex6cap7.txt", h=T)
attach(dados)
dose[1]<- 0.01
ldose<- log10(dose)
 
ajuste1<-glm(as.matrix(dados[,c(2,3)])~ldose, family=binomial(link="logit"), data=dados)
anova(ajuste1,test="Chisq")
ajuste2<-glm(as.matrix(dados[,c(2,3)])~ldose, family=binomial(link="cloglog"), data=dados)
anova(ajuste2,test="Chisq")
ajuste3<-glm(as.matrix(dados[,c(2,3)])~ldose, family=binomial(link="probit"), data=dados)
anova(ajuste3,test="Chisq")
ajuste4<-glm(as.matrix(dados[,c(2,3)])~ldose, family=binomial(link="cauchit"), data=dados)
anova(ajuste4,test="Chisq")
 
cbind(ajuste1$aic, ajuste2$aic, ajuste3$aic, ajuste4$aic)
         [,1]    [,2]     [,3]     [,4]
[1,] 24.64445 27.5475 24.95937 26.14106
 
par(mfrow=c(1,4))
ntot<-matrix(0,6,1)
for (i in 1:6){ntot[i]<-sum(dados[i,2:3])}
 
fit.model<-ajuste1
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste2
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste3
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste4
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
 
## Selecionado o modelo de regressão binomial com link logito (ajuste 1))
## Modelo selecionado (ajuste 1))

ajuste1<-glm(as.matrix(dados[,c(2,3)])~ldose, family=binomial(link="logit"), data=dados)
summary(ajuste1)

dev<-residuals(ajuste1, type='deviance')
dev
QL<-sum(dev^2)
p1<-1-pchisq(QL,4)
cbind(QL,p1)

rpears<-residuals(ajuste1, type='pearson')
rpears

QP<-sum(rpears^2)
p2<-1-pchisq(QP,4)
cbind(QP,p2)

x<-seq(-2,1.5,0.1)
m1<-exp(-4.8869+7.1462*x)/(1+exp(-4.8869+7.1462*x))
par(mfrow=c(1,1))
plot(ldose,sim/(sim+nao),pch=16,ylab="proporção de mortes", xlab="log10(dose)", xlim=c(-2.1,1.5),
     ylim=c(0,1.05))
lines(x,m1, lty=1,lwd=2, col=2)

Com base no modelo ajustado, tem-se, por exemplo, que a chance de morte dos insetos sob a dose de inseticida 5.1 (log10(dose) = 0.7075702) foi exp(7.1462*(0.7075702-0.5797836)) \(\approx\) 2,5 vezes a daqueles sob a dose 3.8 (log10(dose) = 0.5797836).

(b) DL50 = 4.83 e DL90 = 9.8

 10^(-ajuste1$coef[1]/ajuste1$coef[2])
 [1] 4.828918
 
 10^((log(0.9/0.1) + 4.8869)/7.1462)
 [1] 9.801973

Exercício 7

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex7cap7.txt", h=T)
attach(dados)
 
## Obs: Efeito da interação dieta*medicamento não significativo
 
ajuste1<-glm(as.matrix(dados[,c(3,4)])~dieta+medicamento, family=binomial(link="logit"), data=dados)
ajuste2<-glm(as.matrix(dados[,c(3,4)])~dieta+medicamento, family=binomial(link="cloglog"), data=dados)
ajuste3<-glm(as.matrix(dados[,c(3,4)])~dieta+medicamento, family=binomial(link="probit"), data=dados)
ajuste4<-glm(as.matrix(dados[,c(3,4)])~dieta+medicamento, family=binomial(link="cauchit"), data=dados)
cbind(ajuste1$aic, ajuste2$aic, ajuste3$aic, ajuste4$aic)

par(mfrow=c(1,4))
ntot<-matrix(0,9,1)
for (i in 1:9){ntot[i]<-sum(dados[i,3:4])}
fit.model<-ajuste1
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste2
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste3
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
fit.model<-ajuste4
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")

# Obs: devido à facilidade de interpretação optou-se pelo modelo binomial com link logito.
 
anova(ajuste1, test="Chisq")
            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)   
NULL                            8     83.128              
dieta        2   11.125         6     72.004   0.00384 **
medicamento  2   66.743         4      5.261 3.214e-15 ***
dev<-residuals(ajuste1, type='deviance')
QL<-sum(dev^2)
p1<-1-pchisq(QL,4)
cbind(QL,p1)

           QL        p1
[1,] 5.261076 0.2615428
 
summary(ajuste1)
Deviance Residuals:
       1         2         3         4         5         6         7         8         9 
-0.88060  -0.49160   1.36334   0.01514   0.75138  -0.79134   0.85668  -0.19022  -0.65129 
 
Coefficients:
                    Estimate Std. Error z value Pr(>|z|)   
(Intercept)           1.1731     0.1779   6.593 4.30e-11 ***
Dieta_rs             -0.2813     0.1885  -1.492 0.135641   
Dieta_usual          -0.6335     0.1860  -3.405 0.000661 ***
Medicamento_chlor    -0.2724     0.1874  -1.454 0.145958   
Medicamento_placebo  -1.3995     0.1880  -7.445 9.73e-14 ***
 
    Null deviance: 83.1284  on 8  degrees of freedom
Residual deviance:  5.2611  on 4  degrees of freedom
AIC: 58.386

Assim, mantendo-se fixo o medicamento e variando a dieta (sendo a dieta rg = categoria de referência), tem-se:

[i] OR rg | rs = 1/exp(-0.2813) = 1.324851, com IC(OR)_95% = 1/exp(-0.2813 +/- 1.96*0.1885) = (0.916; 1.917)

Embora, nesse estudo, as dietas rs (restrição de sal) e rg (restrição de gordura) não tenham apresentado evidências de diferenças significativas no que se refere à redução da PAD (valor p = 0.135641), notou-se que a chance de redução da PAD sob a dieta com restrição de gordura pode, possivelmente, ser superior a da dieta com restrição de sal, visto que OR rg | rs = 1.32. A realização de estudos adicionais pode auxiliar a avaliar melhor essa questão.

[ii] OR rg | usual = 1/exp(-0.6335) = 1.884194, com IC(OR)_95% = (1.308; 2.713), o que evidencia, para os dados desse estudo, que a chance de redução da PAD sob a dieta com restrição de gordura foi maior do que essa mesma chance sob a dieta usual.

Por outro lado, mantendo-se fixa a dieta e variando o medicamento (atenolol = categoria de referência) tem-se, para esse estudo, que:

[iii] Os medicamentos não apresentaram diferenças estatisticamente significativas (valor p = 0.145958). Contudo, os resultados “sugerem” que a eficácia do chlortalidone, quanto à redução da PAD, talvez seja inferior à do atenolol, uma vez que exp(-0.2724) = 0.7615496. A realização de estudos adicionais pode auxiliar a avaliar melhor essa questão.

[iv] Já OR atenolol | placebo = 1/exp(-1.3995) = 4.053173, com IC(OR)_95% = (2.80; 5.86), o que mostra que a chance de redução da PAD sob o atenolol foi de aproximadamente 4 vezes essa mesma chance sob o placebo.

Capítulo 8 - Regressão multinomial

Exercício 1

[b] Ajustando e selecionando um modelo

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex1cap8.txt", h=T)
attach(dados)
require(VGAM)
 
mlc <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cumulative(parallel=FALSE, reverse=FALSE), dados)
coef(mlc, matrix = TRUE)
 
mop <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cumulative(parallel=TRUE,  reverse=FALSE), dados)
coef(mop, matrix=TRUE)
TRV <- 2*(logLik(mlc)-logLik(mop))
gl <- length(coef(mlc))-length(coef(mop))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
 
         TRV gl           p
[1,] 15.11377  4 0.004470974
 
mlc1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco, cumulative(parallel=FALSE, reverse=FALSE), dados)
mop1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco, cumulative(parallel=TRUE,  reverse=FALSE), dados)
TRV<- 2*(logLik(mlc1)-logLik(mop1))
gl <- length(coef(mlc1))-length(coef(mop1))
p<-1-pchisq(TRV,gl)
cbind(TRV, gl, p)

        TRV gl            p
[1,] 14.88761  2 0.0005850543
 
mlc1 <- vglm(cbind(exc,bom,mod,ruim)~pcard, cumulative(parallel=FALSE, reverse=FALSE), dados)
mop1 <- vglm(cbind(exc,bom,mod,ruim)~pcard, cumulative(parallel=TRUE,  reverse=FALSE), dados)
TRV<- 2*(logLik(mlc1)-logLik(mop1))
gl <- length(coef(mlc1))-length(coef(mop1))
p<-1-pchisq(TRV,gl); cbind(TRV, gl, p)
 
        TRV gl         p
[1,] 0.5879971  2 0.7452776
        
# modelo selecionado 
        
mcpp <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cumulative(parallel=FALSE~tabaco, reverse=FALSE), dados)
coef(mcpp, matrix = TRUE)
 
             logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
(Intercept)     -1.3030321      0.8967713      3.0826522
tabacoSim       -0.1271912     -0.1283847     -0.4581315
pcardSim        -1.0253390     -1.0253390     -1.0253390
 
rp<-resid(mcpp, type = "pearson")
Qp<-sum(rp^2);
cbind(Qp,1-pchisq(Qp,5))

           Qp         
[1,] 6.791648 0.2366026
 
QL<-deviance(mcpp)
cbind(QL,1-pchisq(QL,5))

           QL         
[1,] 6.528623 0.2581221
 
par(mfrow=c(1,2))
rb<- mcpp@y - fitted(mcpp)
plot(1:4,rb[,1], pch=20, ylim=c(-0.5,0.5),xlim=c(1,16),xlab="índice",ylab="Probabilidade observada - predita")
points(5:8, rb[,2], pch=20)
points(9:12, rb[,3], pch=20)
points(13:16, rb[,4], pch=20)
abline(h=0, lty=3)
 
rp<-resid(mcpp, type = "pearson")
plot(1:4,rp[,1], pch=20, ylim=c(-3,3), xlim=c(1,12), xlab="índice", ylab="Resíduos de Pearson")
points(5:8, rp[,2], pch=20)
points(9:12, rp[,3], pch=20)
abline(h=0, lty=3)

[c] Algumas conclusões

A partir do logito 1, tem-se, por exemplo, que a chance de estado de saúde excelente (em relação aos demais) entre os idosos sem problema cardíaco foi 1/exp(-1.025339) = 2.788 vezes a daqueles com problema cardíaco. Já a chance de estado de saúde excelente entre os idosos que não utilizavam tabaco foi 1/exp(-0.1271912) = 1.135 vezes a dos que utilizavam.

Exercício 2

[a] Modelo logitos categorias adjacentes

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/idosos.txt", h=T)
attach(dados)
require(VGAM)
 
mca1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard,acat(parallel=FALSE,reverse=TRUE),dados)
coef(mca1, matrix = TRUE)
 
mca2 <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard,acat(parallel=TRUE,reverse=TRUE),dados)
coef(mca2, matrix=TRUE)
TRV <- 2*(logLik(mca1)-logLik(mca2))
gl <- length(coef(mca1))-length(coef(mca2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

          TRV gl            p
[1,] 22.58513  4 0.0001532303
 
mca1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco,acat(parallel=FALSE,reverse=TRUE),dados)
mca2 <- vglm(cbind(exc,bom,mod,ruim)~tabaco,acat(parallel=TRUE,reverse=TRUE),dados)
TRV <- 2*(logLik(mca1)-logLik(mca2))
gl <- length(coef(mca1))-length(coef(mca2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

          TRV gl           p
[1,] 11.08578  2 0.003915192

mca1 <- vglm(cbind(exc,bom,mod,ruim)~pcard,acat(parallel=FALSE,reverse=TRUE),dados)
mca2 <- vglm(cbind(exc,bom,mod,ruim)~pcard,acat(parallel=TRUE,reverse=TRUE),dados)
TRV <- 2*(logLik(mca1)-logLik(mca2))
gl <- length(coef(mca1))-length(coef(mca2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

          TRV gl           p
[1,] 11.53136  2 0.003133257
 
mca <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard,acat(parallel=FALSE,reverse=TRUE),dados)
summary(mca)
coef(mca, matrix = TRUE)
              loge(P[Y=1]/P[Y=2]) loge(P[Y=2]/P[Y=3]) loge(P[Y=3]/P[Y=4])
(Intercept)         -0.84303147          0.70411995           1.7108833
tabacoSim           -0.08578843         -0.01907128          -0.4264225
pcardSim            -0.61955710         -0.81983308          -0.3516252
 
QL<-deviance(mca)
cbind(QL,1-pchisq(QL,3))
           QL          
[1,] 6.370057 0.09493045    
 
par(mfrow=c(1,2))
rb<- mca@y - fitted(mca)
plot(1:4,rb[,1], pch=20, ylim=c(-0.5,0.5),xlim=c(1,16),xlab="índice",ylab="Probabilidade observada - predita")
points(5:8, rb[,2],pch=20)
points(9:12, rb[,3],pch=20)
points(13:16, rb[,4],pch=20)
abline(h=0, lty=3)
 
rp<-resid(mca, type = "pearson")
plot(1:4,rp[,1], pch=20, ylim=c(-3,3),xlim=c(1,12),xlab="índice",ylab="Resíduos de Pearson")
points(5:8, rp[,2],pch=20)
points(9:12, rp[,3],pch=20)
abline(h=0, lty=3)

A partir do logito 3, tem-se, por exemplo, que a chance de estado de saúde moderado (em relação ao estado ruim) entre os idosos sem problema cardíaco foi 1/exp(-0.3516252) = 1.42 vezes a dos com problema cardíaco. Ainda, a chance de estado de saúde moderado (em relação ao estado ruim) entre os idosos que não utilizavam tabaco foi 1/exp(-0.4264225) = 1.53 vezes a dos que utilizavam.

[b] Modelo logitos razão contínua

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/idosos.txt", h=T)
attach(dados)
require(VGAM)
 
mrc1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cratio(parallel=FALSE, reverse=FALSE), dados)
coef(mrc1, matrix = TRUE)
 
mrc2 <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cratio(parallel=TRUE,  reverse=FALSE), dados)
coef(mrc2, matrix=TRUE)
TRV <- 2*(logLik(mca1)-logLik(mca2))
gl <- length(coef(mrc1))-length(coef(mrc2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

          TRV gl          p
[1,] 11.53136  4 0.02119862
 
mrc1 <- vglm(cbind(exc,bom,mod,ruim)~tabaco, cratio(parallel=FALSE, reverse=FALSE), dados)
mrc2 <- vglm(cbind(exc,bom,mod,ruim)~tabaco, cratio(parallel=TRUE,  reverse=FALSE), dados)
TRV <- 2*(logLik(mrc1)-logLik(mrc2))
gl <- length(coef(mrc1))-length(coef(mrc2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

         TRV gl          p
[1,] 9.08337  2 0.01065544
 
mrc1 <- vglm(cbind(exc,bom,mod,ruim)~pcard, cratio(parallel=FALSE, reverse=FALSE), dados)
mrc2 <- vglm(cbind(exc,bom,mod,ruim)~pcard, cratio(parallel=TRUE,  reverse=FALSE), dados)
TRV <- 2*(logLik(mrc1)-logLik(mrc2))
gl <- length(coef(mrc1))-length(coef(mrc2))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)

         TRV gl            p
[1,] 28.21272  2 7.476286e-07
 
mrc <- vglm(cbind(exc,bom,mod,ruim)~tabaco+pcard, cratio(parallel=FALSE, reverse=FALSE), dados)
summary(mrc)
coef(mrc, matrix = TRUE)
        
             logit(P[Y>1|Y>=1]) logit(P[Y>2|Y>=2]) logit(P[Y>3|Y>=3])
(Intercept)           1.302607         -0.5396741         -1.7103210
tabacoSim             0.124586          0.1052933          0.4252645
pcardSim              1.045950          0.8879319          0.3500455
 
QL<-deviance(mrc)
cbind(QL,1-pchisq(QL,3))

          QL    p-value      
[1,] 6.124047 0.1057286
 
par(mfrow=c(1,2))
rb<- mrc@y - fitted(mrc)
plot(1:4,rb[,1], pch=20, ylim=c(-0.5,0.5),xlim=c(1,16),xlab="índice",ylab="Probabilidade observada - predita")
points(5:8,   rb[,2], pch=20)
points(9:12,  rb[,3], pch=20)
points(13:16, rb[,4], pch=20)
abline(h=0, lty=3)
 
rp<-resid(mrc, type = "pearson")
plot(1:4,rp[,1], pch=20, ylim=c(-3,3),xlim=c(1,12),xlab="índice",ylab="Resíduos de Pearson")
points(5:8, rp[,2],pch=20)
points(9:12, rp[,3],pch=20)
abline(h=0,lty=3)

Exercício 4

[a] Ajustando e selecionando um modelo

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/ex4cap8.txt", h=T)
attach(dados)
require(VGAM)

# modelo inicial (mostrou efeito não significativo da covariável poluição do ar)

mlc <- vglm(cbind(n1,n2,n3,n4)~fumo+trab+ar, cumulative(parallel=FALSE, reverse=FALSE), dados)
summary(mlc)

# avaliando a suposição de chances proporcionais
 
mlc <- vglm(cbind(n1,n2,n3,n4)~fumo+trab, cumulative(parallel=FALSE, reverse=FALSE), dados)
coef(mlc, matrix = TRUE)

            logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
(Intercept)      1.6612419     2.84694620      3.8067147
fumoN            0.4316733     0.03705781      0.9731269
fumoS           -1.3995398    -1.72342341     -1.7783494
trabS           -0.8416610    -0.85802821     -0.9043969
 
mop <- vglm(cbind(n1,n2,n3,n4)~fumo+trab, cumulative(parallel=TRUE, reverse=FALSE), dados)
coef(mop, matrix = TRUE)                  

            logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
(Intercept)      1.7008674      2.5821548      3.5007771
fumoN            0.4013478      0.4013478      0.4013478
fumoS           -1.4505706     -1.4505706     -1.4505706
trabS           -0.8591378     -0.8591378     -0.8591378
 
TRV <- 2*(logLik(mlc)-logLik(mop))
gl <- length(coef(mlc))-length(coef(mop))
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
 
         TRV gl          p
[1,] 12.3538  6 0.05452564
 
mlc1 <- vglm(cbind(n1,n2,n3,n4)~fumo, cumulative(parallel=FALSE, reverse=FALSE), dados)
mop1 <- vglm(cbind(n1,n2,n3,n4)~fumo, cumulative(parallel=TRUE,  reverse=FALSE), dados)
TRV<- 2*(logLik(mlc1)-logLik(mop1))
gl <- length(coef(mlc1))-length(coef(mop1))
p<-1-pchisq(TRV,gl); cbind(TRV, gl, p)
 
          TRV gl          p
[1,] 12.70856  4 0.01279123
 
mlc2 <- vglm(cbind(n1,n2,n3,n4)~trab, cumulative(parallel=FALSE, reverse=FALSE), dados)
mop2 <- vglm(cbind(n1,n2,n3,n4)~trab, cumulative(parallel=TRUE,  reverse=FALSE), dados)
TRV<- 2*(logLik(mlc2)-logLik(mop2))
gl <- length(coef(mlc2))-length(coef(mop2))
p<-1-pchisq(TRV,gl); cbind(TRV, gl, p)
 
           TRV gl         p
[1,] 0.2551593  2 0.8802233

# modelo selecionado: mlc com chances proporcionais parciais            
           
mcpp <- vglm(cbind(n1,n2,n3,n4)~fumo+trab, cumulative(parallel=FALSE~fumo, reverse=FALSE), dados)
coef(mcpp, matrix = TRUE)

            logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3])
(Intercept)      1.6661554     2.84414186      3.7810307
fumoN            0.4306765     0.03814593      0.9705752
fumoS           -1.4014119    -1.72401437     -1.7776227
trabS           -0.8548397    -0.85483967     -0.8548397

[b] Adequação do modelo selecionado: MLC com CPP com as covariáveis fumo e poluição no trabalho

rp<-resid(mcpp, type = "pearson")
Qp<-sum(rp^2);
cbind(Qp,1-pchisq(Qp,26))

          Qp    p-value      
[1,] 13.67476 0.9770414
 
QL<-deviance(mcpp)
cbind(QL,1-pchisq(QL,26))

           QL   p-value       
[1,] 13.75799 0.9760529
 
par(mfrow=c(1,2))
rb<- mcpp@y - fitted(mcpp)
plot(1:12,rb[,1], pch=20, ylim=c(-0.5,0.5), xlim=c(1,48), xlab="índice", ylab="Probabilidade observada - predita")
points(13:24, rb[,2], pch=20)
points(25:36, rb[,3], pch=20)
points(37:48, rb[,4], pch=20)
abline(h=0, lty=3)
 
rp<-resid(mcpp, type = "pearson")
plot(1:12, rp[,1], pch=20, ylim=c(-3,3), xlim=c(1,16), xlab="índice", ylab="Resíduos de Pearson")
points(13:24, rp[,2], pch=20)
points(25:36, rp[,3], pch=20)
abline(h=0, lty=3)

[c] Algumas conclusões

A partir do logito 1 (gravidade I versus II ou III ou IV) tem-se, por exemplo, que:

Poluição no trabalho:   OR S|N = exp(-0.855)   =  0.425    
                        OR N|S = 1/exp(-0.855) =  2.351

Logo, a chance de nenhum sintoma de doença respiratória crônica entre os que trabalhavam em local sem poluição foi 2,35 vezes a dos que trabalhavam em local com poluição.

Quanto ao tabaco, sendo a categoria ex-fumante = nível de referência, tem-se:

                        OR NF|EX = exp(0.43) = 1.537

de modo que a chance de nenhum sintoma de doença respiratória crônica entre os não fumantes (NF) foi em torno de 1,5 vezes a dos ex-fumantes (EX).

Capítulo 9 - Regressão logística condicional

Exercício 1

[a] Análise dos dados

match<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/infart.txt", h=T)
attach(match)
require(survival)
model1<-clogit(MI~SMK+SBP+ECG+SMK*SBP+SMK*ECG+SBP*ECG+strata(par), data=match)
model1
             coef exp(coef) se(coef)     z    p
SMK     -6.81713   0.00109  7.41058 -0.92 0.36
SBP      0.02045   1.02066  0.02740  0.75 0.46
ECG      3.75271  42.63644  7.88380  0.48 0.63
SMK:SBP  0.05828   1.06001  0.05633  1.03 0.30 
SMK:ECG -1.47174   0.22953  1.96530 -0.75 0.45
SBP:ECG -0.00834   0.99169  0.05118 -0.16 0.87
 
Likelihood ratio test=12.7  on 6 df, p=0.0487
n= 78, number of events= 39
 
model2<-clogit(MI~SMK+SBP+ECG+strata(par), data=match)
model2
     coef exp(coef) se(coef)    z     p
SMK 0.837     2.309    0.652 1.28 0.199
SBP 0.035     1.036    0.021 1.66 0.096
ECG 1.831     6.242    1.154 1.59 0.113
 
Likelihood ratio test=11.4  on 3 df, p=0.00963
n= 78, number of events= 39
 
model3<-clogit(MI~SBP+strata(par), data=match)
model3
      coef exp(coef) se(coef)    z     p
SBP 0.0381    1.0388   0.0213 1.79 0.073
ECG 1.6834    5.3837   1.0958 1.54 0.124
 
#  saídas do modelo selecionado

model4<-clogit(MI~SBP+strata(par), data=match)
model4
      coef exp(coef) se(coef)    z     p
SBP 0.0481    1.0492   0.0203 2.36 0.018
 
Likelihood ratio test=6.41  on 1 df, p=0.0113
n= 78, number of events= 39
 
summary(model4)
       coef exp(coef) se(coef)     z Pr(>|z|) 
SBP 0.04807   1.04924  0.02033 2.364   0.0181 *
 
    exp(coef) exp(-coef) lower .95 upper .95
SBP     1.049     0.9531     1.008     1.092
 
plot(model4$residuals, pch=16, ylim=c(-2.5,2.5))

Conclusões com base no modelo selecionado:

  • Para esse estudo, foram encontradas evidências de associação entre a SBP e o infarto do miocárdio (MI), uma vez que a chance de encontrarmos, entre os casos, indivíduos com SBP > 120 mmHg foi maior do que a de encontramos indivíduos com SBP = 120 mmHg.

  • Ou, analogamente, há evidências de associação entre a SBP e o infarto do miocárdio (MI), uma vez que a chance de encontrarmos indivíduos com SBP elevado foi maior entre os casos (indivíduos com MI) do que entre os controles (indivíduos sem MI).

Por exemplo: OR (SBP=140|SBP=120) = exp(0.04807*20) = 2.615. Logo:

  • A chance de encontrarmos, entre os casos, indivíduos com SBP = 140 mmHg foi 2,6 vezes a de encontramos indivíduos com SBP = 120 mmHg. Isso implica que, entre os casos, é esperado 5 indivíduos com SBP = 140 mmHg para cada 2 indivíduos com SBP = 120 mmHg.

  • Ou, analogamente, a chance de encontramos indivíduos com SBP = 140 mmHg entre os casos foi 2,6 vezes a de encontrarmos indivíduos com SBP = 140 mmHg entre os controles. Isso implica que é esperado 5 indivíduos com SBP = 140 mmHg entre os casos para cada 2 com SBP = 140 mmHg entre os controles.

Exercício 3

[a] Analise dos dados

dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/infertilidade.txt", h=T)
i<-order(dados[,7])
dados<-dados[i,]
attach(dados)
 
require(survival)
model1<-clogit(caso~ab_induzido+ab_espontaneo+strata(par), data=dados)
model1
               coef exp(coef) se(coef)    z       p
ab_induzido   1.151     3.162    0.370 3.11  0.0019
ab_espontaneo 1.880     6.551    0.415 4.53 5.9e-06
 
Likelihood ratio test=37.8  on 2 df, p=6.14e-09
n= 166, number of events= 83
 
plot(model1$residuals, pch=149, ylim=c(-2.5,2.5))

[b] Conclusões a partir do modelo ajustado aos dados:

  • A partir do modelo ajustado, há evidências de associação entre a infertilidade secundária feminina e abortos prévios, uma vez que a chance de encontramos, entre os casos, mulheres com 1 ou mais abortos induzidos, bem como mulheres com 1 ou mais abortos espontâneo, foi maior do que a chance de encontramos mulheres com nenhum aborto induzido ou espontâneo.

  • Ou analogamente, há evidências de associação entre a infertilidade secundária feminina e abortos prévios, uma vez que a chance de abortos induzidos, bem como de abortos espontâneos, foi maior entre os casos do que entre os controles.

[1] Variando o número de abortos induzidos (mantendo fixo o número de abortos espontâneos)

  • Entre os casos, notou-se que a chance de encontrarmos mulheres com 1 aborto induzido foi em torno de 3 vezes a de encontramos mulheres com 0 (nenhum) aborto induzido.

  • Ou, analogamente, a chance de mulheres com 1 aborto induzido entre os casos foi em torno de 3 vezes a de mulheres com 1 aborto induzido entre os controles.

[2] Variando o número de abortos espontâneos (mantendo fixo o número de abortos induzidos)

  • Entre os casos, notou-se que a chance de encontramos mulheres com 1 aborto espontâneo foi em torno de 6,5 vezes a de encontrarmos mulheres com 0 (nenhum) aborto espontâneo.

  • Ou, analogamente, a chance de mulheres com 1 aborto espontâneo entre os casos foi em torno de 6,5 vezes a de mulheres com 1 aborto espontâneo entre os controles.