We present a case study with penalized maximum likelihood (PML) estimation of interaction terms. The issue is that assessment of interactions may lead to a large number of formal and informal tests, and iterative adaptations in a prediction model. The model will then likely be overfitted to the data. Shrinkage techniques may reduce the regression coefficients of selected interaction terms such that better predictions are obtained for future patients.
Patients with chronic hepatitis B virus (HBV) infection can be treated with a specific type of interferon (with peginterferon alpha-2b). This leads to response (“HBeAg loss at 26 weeks”) in about a third of the patients. Buster et al aimed to predict response to treatment in individual patients with readily available patient characteristics.
A main effects model was considered with 5 predictors, including a predictor indicating the genotype (A, B, C, D). The researchers were specifically interested in the interaction between genotype (4 groups) and other predictors. When we try perform an overall test for interaction, we encounter a zero cell problem for one interaction term. An overall test for genotype by predictor interaction for 3 of the 4 other predictors is statistically non-significant (p=0.48). But an interaction between genotype and a liver enzyme (“ALT”) is statistically significant if considered in isolation.
The effect of the genotype interaction is that ALT has no effect in genotype A and D, a small effect in genotype C, and a strong effect in genotype B (Fig 12.E1).
Fig 12.E1: Interaction effects in a prediction model for response in HBV patients, with standard and penalized maximum likelihood estimation (see Chapter 13). We note that the ALT effect is especially strong in genotype B patients, but that this effect is attenuated with penalized ML
For estimation, we further consider a penalized maximum likelihood procedure, with separate penalties for main effects and interactions (see Chap 13). This leads to an interesting pattern for the optimal penalty combination. The expected future performance (according to Akaike’s Information Criterion, AIC) is better if interaction terms are penalized more than the main effects (Fig 12.E2).
The optimal combination is a penalty of 0.2 for main effects and 1.2 for interaction terms. Since the interaction was actually selected from a larger number of interaction terms, we might also search for the optimal penalty combination in a model with 3 interactions. This appears to lead to the combination 0, 1.4, implying no penalty for the main effects. An intermediate choice for the model with 1 interaction is to use the penalty combination 0.1, 1.4. This leads to slightly less extreme estimates for the interaction effects, especially for the ALT effect among genotype B patients (Fig 12.E1).
Fig 12.E2: Performance according to AIC for combinations of penalty factors for main effects and an interaction term in the prediction model for response in HBV patients. Better performance is indicated by a larger triangle. The largest triangle is found for the combination 0.2, 1.2 in the model with one interaction term for genotype*ALT
For performance, we consider the c statistics for the main effect model, and for the 2 variants of the interaction model (standard and penalized ML, Table 12.E1. The main effects model has an apparent c statistic of 0.69. With bootstrap validation, the c statistic becomes 0.65. The decrease by 0.04 is substantial, which is explained by the ratio of events per variable (84 events, 5 predictors with 7 df). When we extend the model with the interaction between genotype and ALT, the apparent c increases to 0.72. But the optimism is 0.05, making the model with interaction only slightly better than a main effects model. Using penalized ML led to very similar results, with a slightly lower apparent c statistics and slightly higher validated c statistics. The validation did however not include the search for potential interaction terms, and the search for a reasonable penalty factor. The final model can be presented as a nomogram, which allows for clear insight in the effect of the interaction between genotype and ALT (Fig 12.E3, see Chap 18).
Nomogram presentation of the prediction model for response in HBV patients with penalized estimation of the interaction between ALT and Genotype. For calculation of the probability of response: Determine each predictor value, and draw a straight line upwards to the ‘Points’ axis to determine how many points towards response the patient receives. Sum the points received for each predictor and locate this sum on the ‘Total Points’ axis. Draw a straight line down to find the patient’s predicted probability of response; if the Total Points are less than 3, the probability of response is less than 20%; if over 10, the probability is over 90%.
The essential R code is as follows for penalized estimation of interaction in hepatitis B virus response prediction example.
library(Design) # Fit a model with 3 interaction term in HBVs data set; 5 predictors fulli2 <- lrm(HBEAG78X ~ as.factor(Genotype)*(HBVload+ALT+GGT)+PreviousTx, data=HBVs) # Overall test for interaction anova(fulli2) # Fit a model with 1 interaction term in HBVs data set; 5 predictors fulli <- lrm(HBEAG78X ~ as.factor(Genotype)*ALT + HBVload+GGT+PreviousTx, data=HBVs) # Penalized ML; search penalty factors # 1 penalty; optimal value 0.3 p <- pentrace(fulli, c(0,.2,.25, .3,.35,.4,.6,.8,1.2,2,4,10)) # separate penalty for interaction term; optimal combination 0.2, 1.2 p2 <- pentrace(fulli, list(simple=c(0,.1,.2,.3),interaction=c(0,.4,.8,1,1.2,1.4,1.6))) # plot penalty combinations plot(p2, which='aic.c') # **Fig 12.E2** # separate penalty for 3 interaction terms; optimal combination 0, 1.4 p3 <- pentrace(fulli2, list(simple=c(0,0.01,0.02), interaction=c(0,1,1.2,1.4,1.6,2,4))) # manually change penalty to combination 0.1, 1.4 choice.pen <- p2$penalty choice.pen[1:2] <- c(0.1,1.4) # fit penalized logistic regression model fulli.pen <- update(fulli, penalty=choice.pen) # plot effects of interaction **Fig 12.E1** plot(fulli.pen, ALT=NA, Genotype=NA) # present as nomogram **Fig 12.E3** nomogram(fulli.pen, fun=plogis, lp=F,lp.at=c(-2,0,2,4), ALT=c(2:6), fun.at=seq(.1,.9,by=.1), funlabel="p(response)", vnames="lab", maxscale=5)