· with(predict(quad_bin, dosex_bin, exp=TRUE), {plot(dose, pred, type=“l”, ylim=c(0, 15), ylab=“cardiovascular disease relative risk”, xlab=“alcohol consumption, grams/day”)
· points(dosex_bin$dose, predict(lin_bin, dosex_bin, exp=TRUE)$pred, type=“l”, lty=3, col=“blue”)
The explanation of the graph command is the same as for the linear model.
However, ylim=c(0, 15) means that the y-axis is displayed only between 0 and 15. These values must be adjusted appropriately to represent the graph properly.
The last points command displays the linear model in the same graph.
The black solid line in the middle indicates the risk, and the upper and lower black dotted lines indicate a 95% CI. The blue dotted line at the bottom indicates the risk of the linear model.
Restricted cubic spline model
We will create a cubic spline model among the non-linear models. The method of calculating the regression coefficient is the same as in the abovementioned linear model, but the dose level must be subdivided beforehand.
■ Division of dose level
The dose levels must be subdivided to examine the risk of a dose-response at each level. What should be noted here is that the doses of individual studies are not the same, and the level setting of the dose group can be a problem because the relationship of the risk to the dose-response is non-linear. Therefore, the best approach to the dose group setting in a non-linear relationship that has been developed until now is to subdivide the levels using RCS [
1-
3].
· library(“rms”)
· knots_bin <- quantile(data_bin$dose, c(.05, .35, .65, .95))
In R, RCS can be used through the “rms” package. After specifying data in the quantile function, the 5%, 35%, 65%, and 95% intervals are divided as vectors. The number and range of the dose levels can be adjusted appropriately. When you set four doses, a total of three dose levels are created for knots_bin (knots-1).
■ Calculation of dose level risk (regression coefficient)
· spl_bin <- dosresmeta(formula=logrr ~ rcs(dose, knots_bin), type=type, id=id, se=se, cases=cases, n=n, data=data_bin)
· summary(spl_bin)
The “dosresmeta” function is used. Enter the dependent variable, logrr, in the formula. For independent variables, enter the dose and the previously divided three dose levels after “~” in the RCS function. Then, enter ID, type, se, cases, and n, and set the data_bin with these variables. This model is set as spl_bin.
The goodness of fit of this model is statistically significant (p=0.038).
The estimated regression coefficient for each dose level can be verified, but interpreting the risk for each dose level is meaningless because the dose levels were arbitrarily divided to test the linearity.
■ Testing the linearity of dose level risk (regression coefficient)
The most important part in an analysis that assumes non-linearity is to test whether the slopes of the regression lines for each dose level have statistically significant differences. If they do not show significant differences, it can be determined that there is linearity.
In an analysis that assumes non-linearity, the differences in the slopes of the three dose levels can be tested as follows:
· waldtest(b=coef(spl_bin), Sigma=vcov(spl_bin), Terms=2:3)
The command for testing linearity is “waldtest”. The first of the three dose levels is excluded because it is the total raw data itself, and the joint slope of the second and third dose levels is tested.
The null hypothesis (H0): doses2=doses3=0.
If the p-value is larger than 0.05, the null hypothesis is accepted, and the joint slope is zero, meaning that there is no slope. Because the slopes of the two dose levels have no difference, it is determined that the model has linearity.
If the p-value is smaller than 0.05, the null hypothesis is rejected, and the joint slope is not zero. Thus, because there is a slope or the slopes of the two dose levels have a difference, it is determined that the model has nonlinearity.
As a result of waldtest, the p-value is 0.018. Therefore, this model has non-linearity.
- Cubic spline model graph
The cubic spline model is plotted in
Figure 5.
· xref_bin <- 0
· with(predict(spl_bin, dosex_bin, xref_bin, exp=TRUE),{plot(get(“rcs(dose, knots_bin)dose”), pred, type=“l”, ylim=c(0.4, 10), ylab=“cardiovascular disease relative risk”, xlab=“alcohol consumption, grams/day”, log=“y”, bty=“l”, las=1) matlines (get(“rcs(dose, knots_bin)dose”), cbind(ci.ub, ci.lb), col=1, lty=“dashed”)})
· points(dosex_bin$dose, predict(lin_bin, dosex_bin, xref_bin, exp=TRUE)$pred, type=“l”, lty=3, col=“blue”)
The explanation for the graph command is the same as in the linear model. xref_bin specifies a certain value for the dose reference variable. If you change the value of the reference variable appropriately according to the shape of the created graph, your discriminating ability will improve. Hence, you can set the value to zero at first and change it according to the shape of the graph.
ylim=c(0.4, 10) means that the y-axis is displayed only between 0.4 and 10. The values must be adjusted appropriately to represent the graph properly. The last “points” command displays the linear model in the same graph.
The black solid line in the middle indicates the risk, and the upper and lower black dotted lines indicate a 95% CI. The blue dotted line at the bottom indicates the risk of the linear model.
The risk in the cubic spline model is slightly different from the risk in the linear model. It is similar until a dose of approximately 25, but the risk tends to increase after that and increases sharply from 40. Therefore, according to the cubic spline model, when the alcohol intake dose is high (higher than 30), the risk of CVD increases sharply (
Figure 5).
- Adjusting the cubic spline model graph
In
Figure 5, the width of the 95% CI begins to decrease at a dose of approximately 17 and is recognized as an inflexion point. This value can be checked as follows:
· pre_bin <- predict(spl_bin, dosex_bin, exp=TRUE)
· pre_bin$ci.ub - pre_bin$ci.lb
When you calculate the width of the CI after obtaining all of the predicted values of the RCS model, you can see that the width is the smallest at the 17th dose, which is 0.2679, in the console window. Thus, the graph can be drawn again by setting this as the inflection point, as follows:
· xref_bin <- 17
When the command for creating the cubic spline model graph is performed after changing the dose reference variable to 17, the graph in
Figure 6 is obtained.
After changing the value of the reference variable,
Figure 6 shows a J shape, which rises at the left and right around the dose of 17. Therefore, more intuitive interpretation is possible.