I am fairly new to statistical analysis and was hoping to get some advice on an analysis I am hoping to run.

I have data for children with acute kidney injuries (AKI) classified as a multilevel categorical variable as so:


I also have multiple binary (Y/N or 1/0) variables for these children such as: prematurity, exposure to nephrotox drugs etc.

I was hoping to look for the statistical significance of these factors on AKI. From simple box plots several seem significant. I also ran fisher tests/ chi squared tests on each which show significance in some.

I was hoping to shore-up my analysis and be able to state which variables are independent risk factors but I am nervous I am drawing incorrect conclusions about the data. I would ultimately like to have relative risk scores.

Any information or tests to run would be much appreciated.

Many thanks, Jack


You are likely looking for Ordered Logistic Regression (also called Probit Regression). This allows you to model an ordered response variable (like degree of acute kidney injury) without having to specify the relationship between them (e.g., the step between level 1 and level 2 can be different than the step between level 2 and level 3).

This is relatively easy with polr from the MASS package in R.

First, I am generating some sample data (using tidyverse functions):

n_inds <- 1000

for_factors <-
    , function(idx){
      this_prob <- runif(1)
      sample(0:1, n_inds, replace = TRUE, prob = c(1-this_prob, this_prob))
    }) %>%
  set_names(paste0("factor_", 1:length(.))) %>%

test_data <-
    ind = 1:n_inds
  ) %>%
  ) %>%
    risk_ratio = 
      0.25 +
      factor_1 * .25 +
      factor_2 * .5 +
      factor_3 * .33 +
      factor_4 * .75 +
      factor_5 * .1
    , prob_inj = 1 / ( (1/risk_ratio) + 1)
    , AKI_score = sapply(prob_inj, function(p){
      sum(sample(0:1, 3, replace = TRUE, prob = c(1-p, p)))
    , AKI_score = factor(AKI_score)

My random generation gave me this distribution of scores:

  0   1   2   3 
132 375 345 148 

For simplicities sake, I am creating a second data frame with just the variables that I want to model:

for_model_data <-
  test_data %>%
  dplyr::select(AKI_score, starts_with("factor_"))

Then, I run the model (being able to use . is really the only advantage of making the separate data frame.

aki_model <-
  polr(AKI_score ~ .
       , data = for_model_data
       , Hess = TRUE)

We can evaluate this with summary to show the modeled effects and the thresholds for moving between the categories:

polr(formula = AKI_score ~ ., data = for_model_data, Hess = TRUE)

           Value Std. Error t value
factor_1  0.1561     0.2345  0.6654
factor_2  0.6617     0.1487  4.4501
factor_3  0.5954     0.1201  4.9562
factor_4  1.1735     0.1230  9.5432
factor_5 -0.1307     0.1582 -0.8258

    Value   Std. Error t value
0|1 -1.1301  0.2797    -4.0404
1|2  0.9531  0.2770     3.4414
2|3  2.8462  0.2907     9.7905

Residual Deviance: 2429.711 
AIC: 2445.711 

Note that this does not give significance values. However, we can get confidence intervals from confint(aki_model):

              2.5 %    97.5 %
factor_1 -0.3040116 0.6165130
factor_2  0.3708770 0.9540038
factor_3  0.3604931 0.8315159
factor_4  0.9337777 1.4159578
factor_5 -0.4410188 0.1795969

You will note that factors 1 and 5 (the smallest effects from the model that I built) have confidence intervals that include 0. This is just due to a lack of power to detect such small effects in the model (especially when there are larger effects in there masking the signal).

We can formalize the evaluation of which factors to include by using step-wise model selection with

stepped_model <-

Which returns a model that only includes factors 2, 3, and 4.

This is likely what you want to use if you are not sure which of your factors have an impact. However, note that it excluded two that had a real effect here (we know because we built the model). Further, multiple testing will sometimes fool you. I ran the same model with 50 factors (only the first 5 had a modeled effect) and the stepwise modeling included factors 2, 4, 12, 20, 21, 27, 36, 39, 45, and 49. It excluded three true effects and included eight factors that had no effect due to over-fitting. So, be careful, and take some of these results with a grain of salt.

  • $\begingroup$ Thank you so much this is brilliant! I will update with you with how I get on. $\endgroup$ – Jack Ingham Nov 18 '19 at 14:57
  • $\begingroup$ I'm glad that it helps. If it answers your question (I believe it does), please click check mark (just below the question score) to indicate that you are accepting this answer. It helps other users know that your question has been answered. $\endgroup$ – Mark Peterson Nov 18 '19 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.