Skip to contents

Poisson and Logistic Residuals

Source: vignettes/Residual.Rmd
Residual.Rmd
Sys.setenv(OMP_THREAD_LIMIT = 1) # Reducing core use, to avoid accidental use of too many cores
library(Colossus)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%
if (system.file(package = "survival") != "") {
  library(survival)
}
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:data.table':
#> 
#>     between, first, last
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

Rate and Probability Residuals

After running a regression, someone’s first question might be the accuracy of the results. This vignette will discuss the three general residual options, their equations for each regression type, and how the residual function can be called.

Calculating Residuals for Poisson Models

Each residual (rr) is written in terms of the events (yy), duration (pp), and predicted per duration rate (λ\lambda) for each row (ii). Similar to the relative risk functions, the residual functions were designed to be used with either a training or testing data set. This means that they do not filter the input data, so the residuals are also valid for data that would typically be filtered during a regression. One example would be strata without events or rows without events or duration, which would all be filtered for forcing the predicted event rate to zero.

The first option is the raw residual. This is equal to the difference between the observed and predicted numbers of events. This is the simplest residual, and can be used to visualize simple trends within bins.

ri=yi−λi⋅pi \begin{aligned} r_i = y_i - \lambda_i \cdot p_i \end{aligned}

The second option is the deviance residual. This is equal to the contribution of each row to the deviance. Similar to the standard deviance calculation, the logarithm product is set to zero for rows with no events.

ri=2(yi⋅log(yiλi⋅pi)−(yi−λi⋅pi))ri(yi=0)=0−2(yi−λi⋅pi) \begin{aligned} r_i = 2 \left (y_i \cdot \log \left ( \frac{y_i}{\lambda_i \cdot p_i} \right ) - \left (y_i - \lambda_i \cdot p_i \right ) \right) \\ r_i(y_i = 0) = 0 - 2\left (y_i - \lambda_i \cdot p_i \right) \end{aligned}

The third option is the pearson residual. In the case that both the observed and predicted events are zero, the residual is also set to zero.

ri=(yi−λi⋅pi)2/(λi⋅pi)ri(yi=0,λi⋅pi=0)=0 \begin{aligned} r_i = \left (y_i - \lambda_i \cdot p_i \right )^2 / \left (\lambda_i \cdot p_i \right) \\ r_i(y_i = 0, \lambda_i\cdot p_i = 0) = 0 \end{aligned}

Calculating Residuals for Logistic Regressions

Very similar residual options are available for logistic regressions and models. In this case, the residuals (rr) are written in terms of the observed events (yy), number of trials (nn), and predicted probability (pp). Once again, the first residual is the raw residual and the difference between the observed and predicted events.

ri=yi−pi⋅ni \begin{aligned} r_i = y_i - p_i \cdot n_i \end{aligned}

The second option is the deviance residual, which calculates the contribution of each row to the deviance. Similar to the deviance residual for poisson models, when the number of events are equal to zero or equal to the number of trials, the corresponding logarithm products are set to zero.

ri=|yi−pi⋅ni|yi−pi⋅ni2[yilog(yipi⋅ni)+(ni−yi)log(ni−yini−ni⋅pi)]ri(yi=0)=−2[nilog(11−pi)]ri(yi=ni)=2[yilog(1pi)] \begin{aligned} r_i = \frac{\left| y_i - p_i \cdot n_i \right|}{y_i - p_i \cdot n_i} \sqrt{2 \left[ y_i \log \left ( \frac{y_i}{p_i \cdot n_i} \right ) + \left ( n_i - y_i \right) \log \left ( \frac{n_i - y_i}{n_i - n_i \cdot p_i} \right) \right]}\\ r_i(y_i = 0) = -\sqrt{2 \left[ n_i \log \left ( \frac{1}{1 - p_i} \right) \right]} \\ r_i(y_i = n_i) = \sqrt{2 \left[ y_i \log \left ( \frac{1}{p_i} \right ) \right]} \end{aligned}

The final option is the pearson residual. Once again, if both the observed and predicted events are zero, the residual is also set to zero.

ri=yi−pi⋅nini⋅pi⋅(1−pi)ri(yi=0,pi⋅ni=0)=0 \begin{aligned} r_i = \frac{y_i - p_i \cdot n_i} {\sqrt{ n_i \cdot p_i \cdot \left(1 - p_i \right) }} \\ r_i(y_i = 0, p_i \cdot n_i = 0) = 0 \end{aligned}

Running the functions

Residuals for poisson and logistic models and regressions can be calculated using the Residual() function. Note that this is currently only used for poisson and logistic models, Cox residuals can be calculated using the plotSchoenfeld() and plotMartingale() functions. This function can be applied to the results of a regression, or a model object with a parameter guess.

if (system.file(package = "survival") != "") {
  data(reliability, package = "survival")
  capacitor |> setDT()
  df <- copy(capacitor)
} else {
  voltage <- c(200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350)
  temperature <- c(170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180)
  time <- c(439, 904, 1092, 1105, 572, 690, 904, 1090, 315, 315, 439, 628, 258, 258, 347, 588, 959, 1065, 1065, 1087, 216, 315, 455, 473, 241, 315, 332, 380, 241, 241, 435, 455, 1105, 1105, 1105, 1105, 1090, 1090, 1090, 1090, 628, 628, 628, 628, 588, 588, 588, 588, 1087, 1087, 1087, 1087, 473, 473, 473, 473, 380, 380, 380, 380, 455, 455, 455, 455)
  status <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
  df <- data.table(
    voltage = voltage,
    temperature = temperature,
    time = time,
    status = status
  )
}

df$voltage <- (df$voltage - 200) / 150
df$temperature <- (df$temperature - 170) / 10
df$time <- (df$time - 216) / (1105 - 216)

control <- list(ncores = 1, maxiter = 100, verbose = 2)

a_n <- c(-2, 0.01, 0.01)

poismodel <- get_form(Pois(time, status) ~ loglinear(CONST, temperature, voltage), df)$model

poisres <- PoisRun(poismodel, df,
  a_n = a_n, control = control
)

print(poisres)
#> |--------------------------------------------------------------------------------|
#> Final Results
#>      Covariate Subterm Central Estimate Standard Error 95% Confidence Interval
#>         <char>  <char>            <num>          <num>                  <char>
#> 1:       CONST  loglin           -0.779          0.344       (-1.454 - -0.105)
#> 2: temperature  loglin            0.458          0.360        (-0.247 - 1.163)
#> 3:     voltage  loglin            1.442          0.451         (0.558 - 2.325)
#>    2-tail p-value
#>             <num>
#> 1:        0.02355
#> 2:        0.20293
#> 3:        0.00138
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> 
#> Poisson Model Used
#> Person-year Column: 'time'
#> Event Column: 'status'
#> |- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -|
#> -2*Log-Likelihood: 137.99,  Deviance: 75.99,  AIC: 81.99,  BIC: 150.419
#> Iterations run: 5
#> maximum step size: 1.109e-03, maximum first derivative: 5.291e-03
#> Last iteration improved the log-likelihood by: 1.582e-05
#> Analysis converged
#> Records Used: 63, Records Removed: 1
#> Run finished in 0.014 seconds
#> |--------------------------------------------------------------------------------|

There are three options for running the residual function, which are controlled using pearson and deviance boolean input parameters. In any case the rate per duration and raw residuals are calculated and returned, and each parameter can be set to true to calculate the corresponding residual.

res_p <- Residual(poisres, df, pearson = TRUE, deviance = TRUE)
print(res_p$Risk[1:5])
#> [1] 0.4588005 0.4588005 0.4588005 0.4588005 0.7418893
print(res_p$Raw_Residual[1:5])
#> [1] 0.8849128 0.6449328 0.5479087 0.5411995 0.7029105
print(res_p$Pearson_Residual[1:5])
#> [1] 6.8041523 1.1714357 0.6640337 0.6383972 1.6630784
print(res_p$Deviance_Residual[1:5])
#> [1] 2.5543051 0.7810309 0.4919248 0.4758807 1.0216226

a_n <- poisres$beta_0
res_d <- Residual(poismodel, df, a_n = a_n, pearson = TRUE, deviance = TRUE)
print(res_p$Risk[1:5])
#> [1] 0.4588005 0.4588005 0.4588005 0.4588005 0.7418893
print(res_p$Raw_Residual[1:5])
#> [1] 0.8849128 0.6449328 0.5479087 0.5411995 0.7029105
print(res_p$Pearson_Residual[1:5])
#> [1] 6.8041523 1.1714357 0.6640337 0.6383972 1.6630784
print(res_p$Deviance_Residual[1:5])
#> [1] 2.5543051 0.7810309 0.4919248 0.4758807 1.0216226