Beetles

Carregamento de pacotes:

Neste exemplo são utilizados os seguintes pacotes: dplyr, magrittr, tibble, broom, reshape2, ggplot2, plotly, lmtest, ResouceSelection e modEvA. Para a confecção do documento, é necessário o pacote rmarkdown.

# dplyr -- para facilitar manipulacao de dados
if (!require('dplyr')) install.packages('dplyr'); library('dplyr')

# magrittr -- para facilitar manipulacao de dados
if (!require('magrittr')) install.packages('magrittr'); library('magrittr')

# tibble -- para facilitar manipulacao de dados
if (!require('tibble')) install.packages('tibble'); library('tibble')

# reshape2 -- para facilitar manipulacao de dados
if (!require('reshape2')) install.packages('reshape2'); library('reshape2')

# broom -- para facilitar manipulacao de dados
if (!require('broom')) install.packages('broom'); library('broom')

# ggplot2 -- para os plots
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')

# plotly -- para os plots
if (!require('plotly')) install.packages('plotly'); library('plotly')

# lmtest -- para os testes de razao de verossimilhancas e para o teste de wald
if (!require('lmtest')) install.packages('lmtest'); library('lmtest')

# ResourceSelection -- para utilizar uma versão do teste de Hosmer Lemeshow
if (!require('ResourceSelection')) install.packages('ResourceSelection'); library('ResourceSelection')

# modEVA -- para compelemntar a funcao alternativa do teste de Hosmer Lemeshow
if (!require('modEvA')) install.packages('modEvA'); library('modEvA')

Carregamento dos dados:

dose <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)

killed <- c(6, 13, 18, 28, 52, 53, 61, 59)

total <- c(59, 60, 62, 56, 63, 59, 62, 60)

# criando uma coluna do outro resultado: x = 1 (survived); x = 0 (killed)
beetles <- tibble::data_frame(dose, killed, total) %>%
  dplyr::mutate(survived = total - killed, prop = killed/total)

beetles

Plot:

plotly_pallette <- c('#1F77B4', '#FF7F0E', '#2CA02C', '#D62728')

p1 <- beetles %>%
  ggplot()+
  geom_point(aes(x=dose, y=prop), colour = plotly_pallette[1])+
  labs(x="Dose", y="Prop")+
  ggtitle("Beetle Mortality Data")+
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  guides(color=guide_legend(title=NULL))+
  theme(plot.title = element_text(hjust=0.5))

plotly::ggplotly(p1)

Ajuste do modelo:

resp <- as.matrix(
  beetles %>%
  dplyr::select(killed, survived)
)

fit <- glm(resp ~ dose, family=binomial(link="logit"))

summary(fit)

## 
## Call:
## glm(formula = resp ~ dose, family = binomial(link = "logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5139  -0.2808   0.5536   1.0175   1.3788  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -59.188      5.053  -11.71   <2e-16 ***
## dose          33.401      2.839   11.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 274.8744  on 7  degrees of freedom
## Residual deviance:   8.6398  on 6  degrees of freedom
## AIC: 40.821
## 
## Number of Fisher Scoring iterations: 4

-2logLik:

# -2logLik
-2*logLik(fit)

## 'log Lik.' 36.82103 (df=2)

Int. Conf. OR:

# OR
exp(confint.default(fit)) # o metodo *.default* assume normalidade assintotica, por isso eh utilizado.

##                    2.5 %       97.5 %
## (Intercept) 9.864035e-31 3.947270e-22
## dose        1.226916e+12 8.368858e+16

Teste de Razão de Verossimilhanças:

# 1) razao de verossimilhancas
#library(lmtest)
lmtest::lrtest(fit)

Teste de Wald:

# 2) teste de Wald
lmtest::waldtest(fit)

Teste Score:

# 3) Score
anova(fit,test="Rao")

Teste de Hosmer Lemeshow:

No R, há métodos diferentes que implementam o teste de Hosmer Lemeshow. Abaixo, são apresentadas duas formas.

A primeira delas seria utilizando o modelo a partir do qual os dados foram ajustados na forma de tabela de frequências. Mas como nesse caso só há 8 observações, deve-se ajustar o parâmetro g (número de classes/intervalos) da função para um valor menor.

ResourceSelection::hoslem.test(fit$y, fit$fitted.values, g=5)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  fit$y, fit$fitted.values
## X-squared = 0.092295, df = 3, p-value = 0.9927

A segunda forma é ajustando um novo modelo com base nos dados expandidos com a variável resposta na forma de zeros e uns (0,1). Assim ganham-se graus de liberdade e, consequentemente mais “espaço” para os intervalos. Note que na saída do modelo, tem-se exatamente os mesmos valores de coeficientes e de p-valor, mas valores diferentes para AIC, deviance etc. Essa diferença deve-se justamente ao maior número de observações que se tem agora. Na verdade, com essa configuração de dados, todas as sapidas (AIC, etc) passam a bater exatamente com as saídas do SAS.

beetles2 <- as.data.frame(beetles)

# antes do teste, eh necessario transformar para zero e uns
# o nao evento eh sobreviver
response_l <- lapply(1:nrow(beetles2), function(i){
  killed <- rep(c(0,1), c(beetles2[i,"survived"], beetles2[i,"killed"])) # response
  dose <- rep(beetles2[i,"dose"], length(killed)) # variable
  cbind(dose, killed)
})

beetles_newdata <- as.data.frame(Reduce(rbind, response_l))

fit2 <- glm(killed ~ dose, data=beetles_newdata, family=binomial(link="logit"))

# os deviance, df e AIC diferentes enm relacao ao modelo 1 
#...devem-se ao nro de linhas (obs) no modelo expandido
summary(fit2) 

## 
## Call:
## glm(formula = killed ~ dose, family = binomial(link = "logit"), 
##     data = beetles_newdata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7421  -0.6070   0.2171   0.4658   2.3584  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -59.188      5.053  -11.71   <2e-16 ***
## dose          33.401      2.839   11.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 646.28  on 480  degrees of freedom
## Residual deviance: 380.05  on 479  degrees of freedom
## AIC: 384.05
## 
## Number of Fisher Scoring iterations: 5

Para o teste, foi feita uma pequena aleração na função HLfit() do pacote modEvA. O objetivo é replicar os resultados apresentados no output do SAS para o teste de Hosmer Lemeshow. Os ajustes são para solucionar um bug da função e para fazer com que ela retorne também uma tabela das partições que foram geradas para o teste. Esta função é a única que tem retorno parecido com o do SAS porque permite escolher métodos diferentes para gerar os intervalos nos dados. No caso utilizamos bin.method="quantile".

# library(modEvA)
# alteração da função HLfit:
HLfit2 <- function (model = NULL, obs = NULL, pred = NULL, bin.method, n.bins = 10, fixed.bin.size = FALSE, min.bin.size = 15, min.prob.interval = 0.1, simplif = FALSE) {
  # version 1.5 (24 Jun 2015)
  
  if (!is.null(model)) {
    if (!is.null(obs)) message("Argument 'obs' ignored in favour of 'model'.")
    if (!is.null(pred)) message("Argument 'pred' ignored in favour of 'model'.")
    obs <- model$y
    pred <- model$fitted.values
  }
  
  stopifnot(
    length(obs) == length(pred),
    obs %in% c(0, 1),
    pred >= 0,
    pred <= 1
  )
  
  bins <- getBins(obs = obs, pred = pred, bin.method = bin.method, n.bins = n.bins, fixed.bin.size = fixed.bin.size, min.bin.size = min.bin.size, min.prob.interval = min.prob.interval)
  n.bins <- nrow(bins$bins.table)
  
  # next 4 lines: adapted from hosmerlem function in http://www.stat.sc.edu/~hitchcock/diseaseoutbreakRexample704.txt
  observed <- xtabs(cbind(1 - obs, obs) ~ bins$prob.bin)
  expected <- xtabs(cbind(1 - pred, pred) ~ bins$prob.bin)
  chi.sq <- sum((observed - expected) ^ 2 / expected)
  p.value <- 1 - pchisq(chi.sq, df = n.bins - 2)
  rmse <- sqrt(mean((observed - expected) ^ 2))
  
  if (simplif) return(list(chi.sq = chi.sq, p.value = p.value, RMSE = rmse))
  
  # retiramos plot -- nosso ajuste!
  
  # nosso ajuste: partition table
  partition <- data.frame(cbind(matrix(observed[,2]), matrix(expected[,2])), 
                          cbind(matrix(observed[,1]), matrix(expected[,1])))
  
  colnames(partition) <- c("Observed_Event", "Expected_Event", "Observed_Nonevent", "Expected_Nonevent")
  
  return(list(partition = partition, stats = list(chi.sq = chi.sq, DF = n.bins - 2, p.value = p.value), RMSE = rmse))
}

Resultados do teste de Hosmer-Lemeshow:

out_HL <- HLfit2(model=fit2, bin.method = "quantiles", n.bins = 10) # ok!

tibble::as.tibble(out_HL$partition)

out_HL$stats

## $chi.sq
## [1] 8.433355
## 
## $DF
## [1] 6
## 
## $p.value
## [1] 0.2080418

Ajustes pré-plotagem:

# criando obs artificiais para suavizar a curva
x_limits <- beetles %>%
  dplyr::summarise(min = min(dose), max = max(dose)) # obtendo os limites do eixo x

x <- seq(x_limits[[1]], x_limits[[2]], by=0.001) # criando varios ptos entre esses limites -- by precisa ser bem pequeno aqui

# fazendo a predicao das obs artificiais
# ver sol: https://stackoverflow.com/questions/26694931/how-to-plot-logit-and-probit-in-ggplot2
temp.data = data.frame(dose = x) # se nao passar a coluna com o mesmo nome da variavel nao funciona


# Predict the fitted values given the model and hypothetical data
predicted.data <- as.data.frame(predict(fit, newdata = temp.data, 
                                        type="link", se=TRUE))

# Combine the hypothetical data and predicted values
new.data <- cbind(temp.data, predicted.data)

# Calculate confidence intervals
std <- qnorm(0.95 / 2 + 0.5)
new.data$ymin <- fit$family$linkinv(new.data$fit - std * new.data$se)
new.data$ymax <- fit$family$linkinv(new.data$fit + std * new.data$se)
new.data$fit <- fit$family$linkinv(new.data$fit)  # Rescale to 0-1

Plot:

plotly_pallette <- c('#1F77B4', '#FF7F0E', '#2CA02C', '#D62728')

p2 <- ggplot(beetles, aes(x=dose, y=prop)) 
p2 <- p2 + geom_point(colour = plotly_pallette[1]) + 
  geom_ribbon(data=new.data, 
              aes(y=fit, 
                  ymin=ymin, ymax=ymax), alpha = 0.5, 
              fill = '#FFF0F5')+
  geom_line(data=new.data, aes(y=fit), colour = plotly_pallette[2]) + 
  labs(x="Dose", y="Prop")+
  ggtitle("Beetle Mortality Data")+
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  theme(plot.title = element_text(hjust=0.5))
  

plotly::ggplotly(p2)

model.data <- broom::augment(fit)

model.data <- cbind(model.data, .pearson.resid = fit$residuals)

regMat <- model.data %>%
    dplyr::select(resp, .fitted, .dev.resid = .resid, .hat, .cooksd, .std.resid, .pearson.resid) #%>%

p3 <- regMat %>%
  dplyr::select(.pearson.resid, .dev.resid, .fitted) %>%
  reshape2::melt(id.vars=".fitted") %>%
    ggplot()+
    geom_point(aes(x = .fitted, y = value, colour = factor(variable)), alpha = 0.7)+
    geom_hline(yintercept = 0, alpha = 0.5, size=0.4)+
    ggtitle("Residuals")+
    theme_bw()+
    theme(panel.border = element_blank())+ # para ficar igual o plotly
    guides(color=guide_legend(title=NULL))+
    scale_color_manual(values=plotly_pallette[1:2])+ # NAO FUNCIONA NO GGPLOTLY
    theme(plot.title = element_text(hjust=0.5))

plotly::ggplotly(p3)

# library(plotly)
# library(broom)


p4 <- regMat %>% 
  ggplot()+
  geom_point(aes(x=.fitted, y=.dev.resid), colour = plotly_pallette[1]) + 
  geom_smooth(aes(x=.fitted, y=.dev.resid), method = "loess", size = 0.7, fill = '#FFF0F5')+ 
  labs(x="fitted", y="residuals")+
  # guides(color=guide_legend(title=NULL))+
  ggtitle("Deviance Residuals")+
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  theme(plot.title = element_text(hjust=0.5))


p5 <- regMat %>% 
  ggplot()+
  geom_point(aes(x=.fitted, y=.pearson.resid), colour = plotly_pallette[3]) + 
  geom_smooth(aes(x=.fitted, y=.pearson.resid), method = "loess", size = 0.7, fill = '#FFF0F5')+ 
  labs(x="fitted", y="residuals")+
  # guides(color=guide_legend(title=NULL))+
  ggtitle("Pearson Residuals")+
  theme_bw()+
  theme(panel.border = element_blank())+ # para ficar igual o plotly
  theme(plot.title = element_text(hjust=0.5))

plotly::ggplotly(p4)

plotly::ggplotly(p5)

medinflu <-influence.measures(fit)
inf_measures <- tibble::as.tibble(medinflu$infmat)
inf_measures

inf_measures2 <- inf_measures %>%
  reshape2::melt() #%>% # passando para o formato long para fazermos um facet_wrap

inf_measures2$ind <- 1:nrow(inf_measures)

# primeiro cook.d e leverage

p6 <- inf_measures2 %>%
  dplyr::filter(!stringr::str_detect(variable, "dfb.1_|cov.r")) %>%
  ggplot() +
  geom_point(aes(x = ind, y = value, colour = factor(variable)), size=0.9, alpha=0.7)+
  geom_hline(yintercept = 0, size=0.4, alpha = 0.5)+
  facet_wrap(~variable)+
  guides(color=guide_legend(title=NULL))+
  scale_color_manual(values=plotly_pallette)+
  theme(plot.title = element_text(hjust=0.5))+
  labs(title="Influence Diagnostics", x="index", y="value")
  

plotly::ggplotly(p6)