The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
Nomograms
Nomograms are a graphical representation of a model. For other than the simplest of models, nomograms can better explain what the model is doing, compared to just interpreting coefficients. Each predictor variable (and any interactions) is represented by a horizontal scale.
To use a nomogram to make predictions, you look up the value of the predictor on the scale, then go straight up to the top scale (called points) to identify the number of points. You do this for each of the predictors, then sum the points and look up that value on a scale at the bottom of the plot and use that to predict the response value.
With modern computers we no longer need nomograms to help with the calculations, but they are helpful for opening the black box of the model. Things like transformations, polynomials, or splines of a predictor variable can more easily be interpreted from the plot than from a table of coefficients. Link functions on the response can also be visualized from the plot. Read on for an example to see how these plots can help with interpretation and explanation of models.
This vignette covers the basics and some of the options for creating nomograms using the R6Nomogram package. The idea is to create an initial nomogram, then, since it is an R6 object, use the methods to try different changes and improvements until you produce a good looking nomogram (the first one will probably need improvement).
The automated step work well for models with a predict method that has a terms option, such as lm and glm. Another vignette will cover what to do if there is no predict method, or it does not predict terms.
The Model
We first need a model to visualize. We will use the mtcars data set that comes with R and fit a model to predict mpg (miles per gallon).
A few of the potential predictors make more sense as categorical, rather than numeric, predictors, so we will make a copy of the data frame and change the variables to factors.
mtcars2 <- mtcars
mtcars2$cyl <- factor(mtcars2$cyl)
mtcars2$gear <- factor(mtcars2$gear)
mtcars2$vs <- factor(mtcars2$vs, levels=0:1, labels=c("V", "S"))
fit <- glm(mpg ~ poly(wt, 2) + poly(disp,2) + cyl*gear + vs,
data=mtcars2, family=gaussian(link="inverse"))
summary(fit)
Call:
glm(formula = mpg ~ poly(wt, 2) + poly(disp, 2) + cyl * gear +
vs, family = gaussian(link = "inverse"), data = mtcars2)
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.062011 0.007766 7.985 1.72e-07 ***
poly(wt, 2)1 0.045970 0.022222 2.069 0.0525 .
poly(wt, 2)2 0.026467 0.015388 1.720 0.1017
poly(disp, 2)1 0.039702 0.036696 1.082 0.2928
poly(disp, 2)2 -0.032707 0.017188 -1.903 0.0723 .
cyl6 -0.014901 0.008848 -1.684 0.1085
cyl8 -0.009747 0.011208 -0.870 0.3954
gear4 -0.005664 0.005466 -1.036 0.3132
gear5 -0.010201 0.006419 -1.589 0.1285
vsS -0.001406 0.003784 -0.372 0.7142
cyl6:gear4 0.014192 0.008044 1.764 0.0937 .
cyl8:gear4 NA NA NA NA
cyl6:gear5 0.022982 0.011309 2.032 0.0564 .
cyl8:gear5 0.015078 0.010113 1.491 0.1524
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 5.437256)
Null deviance: 1126.05 on 31 degrees of freedom
Residual deviance: 103.31 on 19 degrees of freedom
AIC: 156.32
Number of Fisher Scoring iterations: 6
I do not claim that this is the best model (or even a good model) for predicting mpg, but it illustrates most of the features of R6Nomogram.
Think about trying to explain the model based on the above summary to a non-statistician.
Using glm with the gaussian family and the “inverse” link means that the model actually fits the reciprocal of mpg, gallons per mile instead of miles per gallon.
First Plot
Since glm models are pretty straight forward, the automated parts of R6Nomogram will work directly, we just create a new object based on the model:
library(R6Nomogram)
n1 <- R6Nomogram$new(fit)
1. Loading newdata
2. Populating Variables
3. Creating x varibles
4. Shifting x variables
5. Scaling x variables
6, Prettying x
Pretyying y
Finishing
The steps taken are printed out. Since it finished we do not have to worry about anything. If there had been something tricky, then we would have seen which step the problem occurred on, rerun everything up to that step, and then manually did the rest.
Now we can call the plot method to see the first iteration of our plot (warning: this will not be pretty, we will improve it).
This gives the basic plot, we can start to see some of the patterns, but there are some labels plotted on top of each other and there are some big gaps. We can improve the plot.
Total Points and Response
The default plot does not show the linear predictor. Sometimes it is distracting, but other times it can grant insight. In our model the linear predictor is the gallons per mile variable and can be interesting. We will tell it to plot the linear predictor and also change the label for it and the response. We need to set the field v.pos.r to NULL so that the correct positions are computed the next time we plot.
n1$plot.lp <- TRUE
n1$v.pos.r <- NULL
n1$lp.lab <- "G/M"
n1$resp.lab <- "M/G"
n1$plot()
Categorical Labels
Next we can see that there are several labels plotting on top of each other on the line for the interaction. This is because “4” is the baseline for cyl and “3” is the baseline for gear, so all the interactions including either of those all have the same value. We will replace one of the values with “Other” and the remaining ones with “” so that they will be blank.
dput(n1$x.pretty.vals$`cyl:gear`)
c("4:3", "6:3", "8:3", "4:4", "6:4", "4:5", "6:5", "8:5")
n1$x.pretty.vals$`cyl:gear` <-
c("Other", "", "", "", "6:4", "", "6:5", "8:5")
n1$plot()
Tick Values
The values where to put the ticks was determined by the pretty and related functions, but they did not take into account the non-linearity. Some of the scales have large blank areas where it would make sense to put more tick marks and labels. Let’s take care of that now.
n1$pretty.y(seq(60, 200, by=10))
n1$pretty('wt', c(1.5, 2, 3, 3.5, 4, 4.25, 4.5,
4.75, 5, 5.25, 5.4))
n1$pretty('disp', c(75, 100, 125, 150, 175, 200, 250,
300, 400, 450))
n1$plot()
Additional Options
There is still overlap on the gallons per mile line, lets set some options to make these look better and also move all of the labels further from the lines and make the tick marks a little longer.
n1$options$text.par[['linear predictor']] <- list(cex=0.7)
n1$options$tik.len <- 0.4
n1$options$txt.pos <- 1.2
n1$options$signif.digits <- 3
n1$plot()
Now it looks much better than when we started.
Predictions
We can do one more thing, demonstrate the predictive idea by giving the plot method a row from our data (or a new, matching data frame).
oldpar <- par(oma=c(0,0,1,0))
n1$plot(predict = mtcars2[1,])
cyl:gear : 39.39687
vs : 3.904236
gear : 12.59639
cyl : 0
disp : 45.10572
wt : 1.293572
Total Points: 102.2968
mtext(rownames(mtcars2)[1], line=1)
n1$plot(predict = mtcars2[20,])
cyl:gear : 0
vs : 0
gear : 12.59639
cyl : 41.36501
disp : 0
wt : 1.730594
Total Points: 55.69199
mtext(rownames(mtcars2)[20], line=1)
n1$plot(predict = mtcars2[15,])
cyl:gear : 0
vs : 3.904236
gear : 28.31842
cyl : 14.30722
disp : 47.63867
wt : 89.31539
Total Points: 183.4839
mtext(rownames(mtcars2)[15], line=1)
Tables
We can approximate the number of points from the graph and see a specific set of values when predicting. But, we may want to have tables with more values that we can look up.
Since we will probably want more dense values for the tables, but we do not want to mess up our plot, we will first clone the object, then redo the “tick” locations, then use the table method to show the tables. These tables can also be saved as a list so that they can be copied to a database or spreadsheet to use in other look-up tools.
n2 <- n1$clone()
n2$pretty.y(seq(60, 200, by=5))
n2$pretty('wt', seq(1.5, 5.4, by=0.1))
n2$pretty('disp', seq(75, 450, by=25))
n2$tables()
$wt
wt Points
1 1.5 5.64567238
2 1.6 4.25317933
3 1.7 3.03549929
4 1.8 2.03406486
5 1.9 1.22862005
6 2.0 0.61905387
7 2.1 0.20997714
8 2.2 0.00000000
9 2.3 -0.01084037
10 2.4 0.17759976
11 2.5 0.56529554
12 2.6 1.15226545
13 2.7 1.93850353
14 2.8 2.92400933
15 2.9 4.10878379
16 3.0 5.49282671
17 3.1 7.07613808
18 3.2 8.85871790
19 3.3 10.84056618
20 3.4 13.02168290
21 3.5 15.40206808
22 3.6 17.98172169
23 3.7 20.76064372
24 3.8 23.73883409
25 3.9 26.91623715
26 4.0 30.29284363
27 4.1 33.86921648
28 4.2 37.64566452
29 4.3 41.62197745
30 4.4 45.79789467
31 4.5 50.17315559
32 4.6 54.74749962
33 4.7 59.52066616
34 4.8 64.49239463
35 4.9 69.66242441
36 5.0 75.03049494
37 5.1 80.59634560
38 5.2 86.35971581
39 5.3 92.32140215
40 5.4 98.49651397
$disp
disp Points
1 75 2.390987
2 100 16.822577
3 125 29.696899
4 150 41.014535
5 175 50.775493
6 200 58.979775
7 225 65.627379
8 250 70.718306
9 275 74.252555
10 300 76.230128
11 325 76.651022
12 350 75.515242
13 375 72.823871
14 400 68.573647
15 425 62.760774
16 450 55.415453
$cyl
cyl Points
1 4 41.36501
2 6 0.00000
3 8 14.30722
$gear
gear Points
1 3 28.31842
2 4 12.59639
3 5 0.00000
$vs
vs Points
1 V 3.904236
2 S 0.000000
$`cyl:gear`
cyl:gear Points
1 Other 0.00000
2 0.00000
3 0.00000
4 0.00000
5 6:4 39.39687
6 0.00000
7 6:5 63.79647
8 8:5 41.85567
$Response
Total Points G/M M/G
1 60 0.0330 30.3
2 65 0.0348 28.7
3 70 0.0366 27.3
4 75 0.0384 26.0
5 80 0.0402 24.8
6 85 0.0420 23.8
7 90 0.0439 22.8
8 95 0.0457 21.9
9 100 0.0475 21.1
10 105 0.0493 20.3
11 110 0.0511 19.6
12 115 0.0529 18.9
13 120 0.0547 18.3
14 125 0.0565 17.7
15 130 0.0583 17.2
16 135 0.0601 16.6
17 140 0.0619 16.2
18 145 0.0637 15.7
19 150 0.0655 15.3
20 155 0.0673 14.9
21 160 0.0691 14.5
22 165 0.0709 14.1
23 170 0.0727 13.8
24 175 0.0745 13.4
25 180 0.0763 13.1
26 185 0.0781 12.8
27 190 0.0799 12.5
28 195 0.0817 12.2
29 200 0.0835 12.0
Model Structure
The nomogram is dependent on how you fit the model. Two different models that make the exact same predictions can result in different nomograms. Changing some options in how you fit the model may make the nomogram easier (or harder) to interpret.
For example, the above model used the default treatment contrasts for cyl and gear with 4 cylinders and 3 gears as the baseline values. Therefore, the interactions involving 4 cylinders or 3 gears all had the same value (which we combined into other above).
We can refit the model with a different set of contrasts and see how that changes the plot.
contrasts(mtcars2$cyl) <- contr.sum(3)
contrasts(mtcars2$gear) <- contr.sum(3)
fit2 <- glm(mpg ~ poly(wt, 2) + poly(disp,2) + cyl*gear + vs,
data=mtcars2, family=gaussian(link="inverse"))
n3 <- R6Nomogram$new(fit2)
1. Loading newdata
2. Populating Variables
Warning: contrasts dropped from factor cyl
Warning: contrasts dropped from factor gear
Warning: contrasts dropped from factor cyl
Warning: contrasts dropped from factor gear
Warning: contrasts dropped from factor cyl
Warning: contrasts dropped from factor gear
Warning: contrasts dropped from factor cyl
Warning: contrasts dropped from factor gear
3. Creating x varibles
4. Shifting x variables
5. Scaling x variables
6, Prettying x
Pretyying y
Finishing
n3$options$tik.len <- 0.4
n3$options$txt.pos <- 1.2
n3$options$signif.digits <- 3
n3$plot.lp <- TRUE
n3$lp.lab <- "G/M"
n3$resp.lab <- "M/G"
n3$options$text.par[['linear predictor']] <- list(cex=0.8)
n3$plot()
Since we have to look at the interaction line in the nomogram anyways, we could further simplify the plot by fitting a combination of cyl and gear instead of separate main effects and interaction.
mtcars2$`cyl:gear` <- interaction(mtcars2$cyl, mtcars2$gear,
sep=':')
fit3 <- glm(mpg ~ poly(wt, 2) + poly(disp, 2) +
`cyl:gear`,
data=mtcars2, family=gaussian(link='inverse'))
n4 <- R6Nomogram$new(fit3, verbose=FALSE)
n4$options$tik.len <- 0.4
n4$options$txt.pos <- 1.2
n4$options$signif.digits <- 3
n4$plot.lp <- TRUE
n4$lp.lab <- "G/M"
n4$resp.lab <- "M/G"
n4$x.y.offsets$`cyl:gear` <- c(-1, 1, 2, -1, 1, -1, 1, -2)
n4$options$text.par[['linear predictor']] <- list(cex=0.8)
n4$options$text.par[['cyl:gear']] <- list(col=c(
c("red", "forestgreen", "blue")[c(1,1,1,2,2,3,3,3)]
))
n4$plot()
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.
