Baseado em:
http://www.ats.ucla.edu/stat/mult_pkg/faq/general/odds_ratio.htm
http://www.ats.ucla.edu/stat/r/dae/logit.htm
## Source: local data frame [6 x 10]
##
## id DG30 cond dose idade npart part_trt score pre cio
## 1 17 vacia cio 600 11 5 41 3.5 0 1
## 2 73 prenhe cio 300 9 6 77 3.0 1 1
## 3 74 vacia cio 300 8 6 55 3.5 0 1
## 4 116 prenhe cio 0 6 3 80 3.5 1 1
## 5 120 prenhe cio 900 6 4 NA 2.5 1 1
## 6 150 prenhe cio 0 8 5 89 3.0 1 1Modelo inicial
##
## Call:
## glm(formula = pre ~ dose, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0718 -0.9986 -0.8609 1.2868 1.5312
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2534919 0.3177783 -0.798 0.425
## dose -0.0006090 0.0006127 -0.994 0.320
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.07 on 94 degrees of freedom
## Residual deviance: 125.07 on 93 degrees of freedom
## AIC: 129.07
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 0.7760860 0.4114627 1.442085
## dose 0.9993912 0.9981650 1.000580
Modelo com interação
Dose * condição(anestro / cio)
- logit(p) = log(p/(1-p))= β0 + β1cio + β2dose + β3ciodose *
## , , cond = anestro
##
## DG30
## dose prenhe vacia
## 0 6 1
## 300 1 3
## 600 0 1
## 900 2 1
##
## , , cond = cio
##
## DG30
## dose prenhe vacia
## 0 10 14
## 300 5 14
## 600 5 12
## 900 7 13
## Source: local data frame [6 x 10]
##
## id DG30 cond dose idade npart part_trt score pre cio
## 1 17 vacia cio 600 11 5 41 3.5 0 1
## 2 73 prenhe cio 300 9 6 77 3.0 1 1
## 3 74 vacia cio 300 8 6 55 3.5 0 1
## 4 116 prenhe cio 0 6 3 80 3.5 1 1
## 5 120 prenhe cio 900 6 4 NA 2.5 1 1
## 6 150 prenhe cio 0 8 5 89 3.0 1 1##
## Call:
## glm(formula = pre ~ cond * dose, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5484 -0.9206 -0.8550 1.4181 1.5386
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.839887 0.735107 1.143 0.2532
## condcio -1.389910 0.821970 -1.691 0.0908 .
## dose -0.001393 0.001557 -0.895 0.3708
## condcio:dose 0.001095 0.001700 0.644 0.5194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.07 on 94 degrees of freedom
## Residual deviance: 121.47 on 91 degrees of freedom
## AIC: 129.47
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 2.3161052 0.58612458 11.624445
## condcio 0.2490978 0.04311703 1.174166
## dose 0.9986080 0.99523642 1.001620
## condcio:dose 1.0010958 0.99780584 1.004714
Para vacas em anestro (cond=0), a ecuação é simples:
logit(p) = log(p/(1-p))= β0 + β2*dose
Para vacas em cio (cond=1), a ecuação é:
logit(p) = log(p/(1-p))= (β0 + β1) + (β2 + β3 )*dose
Então podemos dizer que o coeficiente “dose” é o efeito da dose H quando condição=0 (anestro).
Mais explicitamente, quando o tratamento foi aplicado em vacas em anestro, uma unidade de incremento de dose H implica uma mudança no log odds de -0.0014
Por outro lado, para vacas em cio, um incremento unitario de dose H, significa uma mudança no log odds de (-0.001393 + 0.001095) = -0.000298
Em termos de odds ratios, podemos dizer que o hormonio aplicado em vacas em anestro, o OR é exp(-0.0014) = 0.9986 por um incremento unitario de dose H
OR for vacas em cio is exp(-0.000298) = 0.9997 for a one-unit increase in dose de H.
A razao destes dois OR (cio over anestro) resulta ser o exp(coeficiente) para o termo da interação de dose*condcio: 0.9986/0.9997 = exp(0.0010958) = 1.001096
RR = 0.998608/0.9997
Aplicar o modelo para estimar probabilidades
## Source: local data frame [6 x 10]
##
## id DG30 cond dose idade npart part_trt score pre cio
## 1 17 vacia cio 600 11 5 41 3.5 0 1
## 2 73 prenhe cio 300 9 6 77 3.0 1 1
## 3 74 vacia cio 300 8 6 55 3.5 0 1
## 4 116 prenhe cio 0 6 3 80 3.5 1 1
## 5 120 prenhe cio 900 6 4 NA 2.5 1 1
## 6 150 prenhe cio 0 8 5 89 3.0 1 1## pre
## modpred_pre 0 1
## 0 54 29
## 1 5 7## [1] 0.3578947Dose * Idade
##
## Call:
## glm(formula = pre ~ dose * idade, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4946 -1.0176 -0.7813 1.2771 1.7850
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4832905 0.7817457 -0.618 0.536
## dose -0.0016127 0.0014474 -1.114 0.265
## idade 0.0413736 0.1202360 0.344 0.731
## dose:idade 0.0001648 0.0002140 0.770 0.441
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.07 on 94 degrees of freedom
## Residual deviance: 122.07 on 91 degrees of freedom
## AIC: 130.07
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 0.6167506 0.1279419 2.847598
## dose 0.9983886 0.9954582 1.001184
## idade 1.0422414 0.8205855 1.327654
## dose:idade 1.0001648 0.9997471 1.000601
Dose * número de partos
##
## Call:
## glm(formula = pre ~ dose * npart, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6692 -0.9982 -0.8090 1.2729 1.7855
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0590193 0.5616115 -0.105 0.916
## dose -0.0017032 0.0010600 -1.607 0.108
## npart -0.0541964 0.1389997 -0.390 0.697
## dose:npart 0.0003102 0.0002445 1.269 0.205
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.07 on 94 degrees of freedom
## Residual deviance: 122.17 on 91 degrees of freedom
## AIC: 130.17
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 0.9426886 0.3083562 2.849568
## dose 0.9982983 0.9961361 1.000325
## npart 0.9472461 0.7154677 1.243204
## dose:npart 1.0003102 0.9998435 1.000819
Dose * Dias post parto do suministro do hormonio
##
## Call:
## glm(formula = pre ~ dose * part_trt, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2957 -0.9731 -0.8226 1.2946 1.7855
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.366e+00 1.562e+00 -1.515 0.130
## dose 1.265e-03 2.918e-03 0.433 0.665
## part_trt 2.490e-02 1.812e-02 1.374 0.169
## dose:part_trt -2.046e-05 3.397e-05 -0.602 0.547
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 106.82 on 79 degrees of freedom
## Residual deviance: 103.89 on 76 degrees of freedom
## (15 observations deleted due to missingness)
## AIC: 111.89
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 0.09388739 0.003223408 1.673830
## dose 1.00126544 0.995411924 1.007110
## part_trt 1.02521481 0.991119249 1.065572
## dose:part_trt 0.99997954 0.999911534 1.000047
Dose * Score
##
## Call:
## glm(formula = pre ~ dose * score, family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0769 -0.9852 -0.9026 1.2880 1.6121
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2865823 2.1665995 -0.132 0.895
## dose -0.0016188 0.0043773 -0.370 0.712
## score 0.0121592 0.7111642 0.017 0.986
## dose:score 0.0003254 0.0014138 0.230 0.818
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 126.07 on 94 degrees of freedom
## Residual deviance: 124.94 on 91 degrees of freedom
## AIC: 132.94
##
## Number of Fisher Scoring iterations: 4## OR 2.5 % 97.5 %
## (Intercept) 0.7508253 0.01010977 53.874982
## dose 0.9983825 0.98964591 1.006946
## score 1.0122334 0.24709088 4.134484
## dose:score 1.0003254 0.99754992 1.003147
Modelo completo com efeitos aditivos
Dose como variavel continua
## $first.line
## [1] "\nLogistic regression predicting pre \n"
##
## $table
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test) P(LR-test)
## dose (cont. var.) "0.9996 (0.9983,1.0009) " "0.9996 (0.9981,1.001) " "0.561" "0.56"
## "" "" "" ""
## cond: cio vs anestro "0.33 (0.1,1.11) " "0.24 (0.06,0.93) " "0.038" "0.033"
## "" "" "" ""
## idade (cont. var.) "1.12 (0.97,1.31) " "1.17 (0.99,1.37) " "0.066" "0.064"
## "" "" "" ""
## score (cont. var.) "1.35 (0.5,3.59) " "1.49 (0.51,4.35) " "0.47" "0.469"
## "" "" "" ""
## part_trt (cont. var.) "1.02 (0.99,1.04) " "1.02 (1,1.05) " "0.076" "0.066"
## "" "" "" ""
##
## $last.lines
## [1] "Log-likelihood = -48.237\nNo. of observations = 80\nAIC value = 108.474\n\n"
##
## attr(,"class")
## [1] "display" "list"Dose como fator
##
## Call:
## glm(formula = pre ~ factor(dose) + idade + score + part_trt,
## family = "binomial", data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8457 -0.8968 -0.6483 1.0194 2.4566
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.76831 2.07518 -1.334 0.1822
## factor(dose)300 -1.58694 0.71218 -2.228 0.0259 *
## factor(dose)600 -1.04670 0.69878 -1.498 0.1342
## factor(dose)900 -0.41031 0.64955 -0.632 0.5276
## idade 0.13942 0.08722 1.599 0.1099
## score 0.08871 0.56000 0.158 0.8741
## part_trt 0.02200 0.01298 1.695 0.0901 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 106.819 on 79 degrees of freedom
## Residual deviance: 95.205 on 73 degrees of freedom
## (15 observations deleted due to missingness)
## AIC: 109.21
##
## Number of Fisher Scoring iterations: 4- Por cada cambio unitario em idade, o log odds de ficar prenhe (versus vacia) aumentou em 0.139
- Por cada cambio unitario em score, o log odds de ficar prenhe (versus vacia) aumentou em 0.088
Por cada cambio unitario em dias post parto ao suministro do hormonio, o log odds de ficar prenhe (versus vacia) aumentou em 0.022
- Havendo recebido 300 ml de hormonio (vs nao receber nada), mudou o log odds de “ficar prenhe” por -1.586
- Havendo recebido 600 ml de hormonio (vs nao receber nada), mudou o log odds de “ficar prenhe” por -1.046
Havendo recebido 900 ml de hormonio (vs nao receber nada), mudou o log odds de “ficar prenhe” por -0.41
Testamos a significancia do uso das 3 doses do hormonio.
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 5.8, df = 3, P(> X2) = 0.12O teste estatistico chi-squared 5.8, com 3 graus de liberdade está asociado com o p-valor de 0.12 indicando o efeito geral do uso do hormonio nao foi estatísticamente significativo.
Podemos testar a significancia das doses individuais, ou seja, usar 300, 600 ou 900 ml(?) de hormonio vs nao usar nada,
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 0.27, df = 1, P(> X2) = 0.6## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 0.63, df = 1, P(> X2) = 0.43## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 1.2, df = 1, P(> X2) = 0.28Para interpretar eles como odds-ratio » exp(coeficientes)
## OR 2.5 % 97.5 %
## (Intercept) 0.06276807 0.0009258949 3.4300446
## factor(dose)300 0.20455116 0.0453709952 0.7742001
## factor(dose)600 0.35109529 0.0830598431 1.3312139
## factor(dose)900 0.66344641 0.1802508462 2.3604655
## idade 1.14960691 0.9720288316 1.3751250
## score 1.09276339 0.3580793797 3.2841743
## part_trt 1.02223885 0.9976103943 1.0503940Agora podemos dizer que por um incremento unitario em idade (ano), o odds de ficar prenhe (versus vacia) incrementa por um fator de 1.14 veces…
Probabilidades preditas pelo modelo:
Criamos uma nova tabela de dados variando uma variavel por vez e mantendo as restante nas médias e calculamos a probabilidade predita de ficar prenhe nessas condiçoes (doses de hormônio)
Tabela de probabilidades preditas pelo modelo
In the above output we see that the predicted probability of ficar prenhe is 0.531 for vacas sem uso de hormonio (dose=0), e 0.429 for vacas com maior dose de hormonio (900 ml), holding idade and score at their means.
We can do something very similar to create a table of predicted probabilities varying the value of idade and score.
We are going to plot these, so we will create 100 values of idade between 13 and 14, at each value of dose.
The code to generate the predicted probabilities (the first line below) is the same as before, except we are also going to ask for standard errors so we can plot a confidence interval. We get the estimates on the link scale and back transform both the predicted values and confidence limits into probabilities.
It can also be helpful to use graphs of predicted probabilities to understand and/or present the model. We will use the ggplot2 package for graphing. Below we make a plot with the predicted probabilities, and 95% confidence intervals.
Para score
Para dias post-parto ao suministro dos hormônios
## id DG30 cond dose idade npart
## Min. : 8.0 prenhe:36 anestro:15 Min. : 0.0 Min. : 3.000 Min. : 0.000
## 1st Qu.: 128.0 vacia :59 cio :80 1st Qu.: 0.0 1st Qu.: 3.000 1st Qu.: 1.000
## Median : 181.0 Median :300.0 Median : 5.000 Median : 3.000
## Mean : 278.7 Mean :404.2 Mean : 5.747 Mean : 3.263
## 3rd Qu.: 264.5 3rd Qu.:600.0 3rd Qu.: 7.500 3rd Qu.: 5.000
## Max. :4944.0 Max. :900.0 Max. :14.000 Max. :12.000
##
## part_trt score pre cio
## Min. : 33.00 Min. :2.000 Min. :0.0000 Length:95
## 1st Qu.: 70.25 1st Qu.:2.500 1st Qu.:0.0000 Class :character
## Median : 86.00 Median :3.000 Median :0.0000 Mode :character
## Mean : 82.97 Mean :3.037 Mean :0.3789
## 3rd Qu.: 96.25 3rd Qu.:3.500 3rd Qu.:1.0000
## Max. :120.00 Max. :4.000 Max. :1.0000
## NA's :15