{"id":833,"date":"2023-08-13T18:28:44","date_gmt":"2023-08-13T22:28:44","guid":{"rendered":"https:\/\/www.clayford.net\/statistics\/?p=833"},"modified":"2023-08-19T07:50:29","modified_gmt":"2023-08-19T11:50:29","slug":"rms-for-the-rest-of-us","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/rms-for-the-rest-of-us\/","title":{"rendered":"rms for the rest of us"},"content":{"rendered":"

For a few years now I\u2019ve been chipping away at Regression Modeling Strategies<\/em><\/a> by Frank Harrell<\/a>. Between the exposition, the case studies, the notes, the references, the code and exercises, it\u2019s a deep and rich source of statistical learning. The book\u2019s associated R package, {rms}, is used throughout the book and is one of the reasons it can be tough sledding at times, especially if you\u2019re like me and learned regression using the base R functions lm()<\/code> and glm()<\/code>. It’s not that {rms} is hard to use. In fact it’s quite easy to use. But using functions such as summary()<\/code> and anova()<\/code> with {rms} models produces very different output than what you get with a base R lm()<\/code> model and can seem baffling to the uninitiated. In this post I hope to explain a little of what {rms} is doing by comparing it to the more traditional approaches in R so commonly taught in classrooms and textbooks.<\/p>\n

\n

Model fitting and summary output<\/h2>\n

To begin, let\u2019s load the gala<\/code> data from the {faraway} package. The gala<\/code> data set contains data on species diversity on the Galapagos Islands. Below we use the base R lm()<\/code> function to model the number of plant Species<\/em> found on the island as function of the Area<\/em> of the island (km\\(^2\\)<\/span>), the highest Elevation<\/em> of the island (m), the distance from the Nearest island (km), the distance from Santa Cruz<\/em> island (km), and the area of the Adjacent<\/em> island (square km). This is how I and thousands of others learned to do regression in R. Of course we use the familiar summary()<\/code> function on our model object to see the model coefficients, marginal tests, the residual standard error, R squared, etc.<\/p>\n

library(faraway)\r\ndata("gala")\r\nm <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, \r\n        data = gala)\r\nsummary(m)<\/code><\/pre>\n
## \r\n## Call:\r\n## lm(formula = Species ~ Area + Elevation + Nearest + Scruz + Adjacent, \r\n##     data = gala)\r\n## \r\n## Residuals:\r\n##      Min       1Q   Median       3Q      Max \r\n## -111.679  -34.898   -7.862   33.460  182.584 \r\n## \r\n## Coefficients:\r\n##              Estimate Std. Error t value Pr(>|t|)    \r\n## (Intercept)  7.068221  19.154198   0.369 0.715351    \r\n## Area        -0.023938   0.022422  -1.068 0.296318    \r\n## Elevation    0.319465   0.053663   5.953 3.82e-06 ***\r\n## Nearest      0.009144   1.054136   0.009 0.993151    \r\n## Scruz       -0.240524   0.215402  -1.117 0.275208    \r\n## Adjacent    -0.074805   0.017700  -4.226 0.000297 ***\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\r\n## \r\n## Residual standard error: 60.98 on 24 degrees of freedom\r\n## Multiple R-squared:  0.7658, Adjusted R-squared:  0.7171 \r\n## F-statistic:  15.7 on 5 and 24 DF,  p-value: 6.838e-07<\/code><\/pre>\n
Now let\u2019s fit the same model using the {rms} package. For this we use the ols()<\/code> function. Notice we can use R\u2019s formula syntax as usual. One additional argument we need to use that will come in handy in a few moments is x = TRUE<\/code>. This stores the predictor variables with the model fit as a design matrix. Notice we don\u2019t need to call summary()<\/code> on the saved model object. We simply print it.<\/p>\n
library(rms)\r\nmr <- ols(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, \r\n            data = gala, x = TRUE)\r\nmr<\/code><\/pre>\n
## Linear Regression Model\r\n## \r\n## ols(formula = Species ~ Area + Elevation + Nearest + Scruz + \r\n##     Adjacent, data = gala, x = TRUE)\r\n## \r\n##                  Model Likelihood    Discrimination    \r\n##                        Ratio Test           Indexes    \r\n## Obs       30    LR chi2     43.55    R2       0.766    \r\n## sigma60.9752    d.f.            5    R2 adj   0.717    \r\n## d.f.      24    Pr(> chi2) 0.0000    g      105.768    \r\n## \r\n## Residuals\r\n## \r\n##      Min       1Q   Median       3Q      Max \r\n## -111.679  -34.898   -7.862   33.460  182.584 \r\n## \r\n## \r\n##           Coef    S.E.    t     Pr(>|t|)\r\n## Intercept  7.0682 19.1542  0.37 0.7154  \r\n## Area      -0.0239  0.0224 -1.07 0.2963  \r\n## Elevation  0.3195  0.0537  5.95 <0.0001 \r\n## Nearest    0.0091  1.0541  0.01 0.9932  \r\n## Scruz     -0.2405  0.2154 -1.12 0.2752  \r\n## Adjacent  -0.0748  0.0177 -4.23 0.0003<\/code><\/pre>\n
The coefficient tables and summaries of residuals are the same. Likewise both return R-squared and Adjusted R-squared. The residual standard error from the lm()<\/code> output is called sigma<\/code> in the ols()<\/code> output.<\/p>\n
The ols<\/code> output also includes a \u201cg\u201d statistic. This is Gini\u2019s mean difference and measures the dispersion of predicted values. It\u2019s an alternative to standard deviation. We can use the {Hmisc} function GiniMd<\/code> to calculate Gini\u2019s mean difference for the model fit with lm()<\/code> as follows.<\/p>\n
Hmisc::GiniMd(predict(m))<\/code><\/pre>\n
## [1] 105.7677<\/code><\/pre>\n
Whereas the summary output for lm()<\/code> includes an F test for the null that all of the predictor coefficients equal 0, the ols()<\/code> function reports a Likelihood Ratio Test. This tests the same null hypothesis using a ratio of likelihoods. We can calculate this test for the lm()<\/code> model using the logLik()<\/code> function. First we subtract the log likelihood of a model with no predictors from the log likelihood for the original model, and then multiply by 2. Apparently, according to Wikipedia<\/a>, multiplying by 2 ensures mathematically that the test statistic converges to a chi-square distribution if the null is true. If the null is true, we expect this different of log likelihoods (or ratio, according to the quotient rule of logs<\/a>) to be about 1. This statistic is very large. Clearly at least one of the predictor coefficients is not 0.<\/p>\n
LRchi2 <- (logLik(m) - logLik(lm(gala$Species ~ 1)))*2\r\nLRchi2<\/code><\/pre>\n
## 'log Lik.' 43.55341 (df=7)<\/code><\/pre>\n
If we wanted to calculate the p-value that appears in the ols()<\/code> output, we can use the pchisq()<\/code> function. Obviously the p-value in the ols()<\/code> output is rounded.<\/p>\n
pchisq(LRchi2, df = 7, lower.tail = FALSE)<\/code><\/pre>\n
## 'log Lik.' 2.60758e-07 (df=7)<\/code><\/pre>\n<\/div>\n
\n
Partial F tests and important measures<\/h2>\n
Now let\u2019s use the anova()<\/code> function on both model objects and compare the output. For the lm()<\/code> object, we get sequential partial F tests using Type I sums of squares. Each line tests the null hypothesis that adding the listed predictor to the previous model without it explains no additional variability. So the Area line compares a model with just an intercept to a model with an intercept and Area. The Elevation line compares a model with an intercept and Area to a model with an intercept, Area, and Elevation. And so on. Small p-values are evidence against the null.<\/p>\n
anova(m)<\/code><\/pre>\n
## Analysis of Variance Table\r\n## \r\n## Response: Species\r\n##           Df Sum Sq Mean Sq F value    Pr(>F)    \r\n## Area       1 145470  145470 39.1262 1.826e-06 ***\r\n## Elevation  1  65664   65664 17.6613 0.0003155 ***\r\n## Nearest    1     29      29  0.0079 0.9300674    \r\n## Scruz      1  14280   14280  3.8408 0.0617324 .  \r\n## Adjacent   1  66406   66406 17.8609 0.0002971 ***\r\n## Residuals 24  89231    3718                      \r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1<\/code><\/pre>\n
Calling anova()<\/code> on the ols()<\/code> object produces a series of partial F tests using Type II sums of squares. Each line tests the null hypothesis that adding the listed predictor to a model with all the other listed predictors already in it explains no additional variability. So the Area line compares a model with all the other predictors to a model with all the other predictors and Area. The Elevation line compares a model with all the other predictors to a model with all the other predictors and Elevation. The line labeled REGRESSION is the F Test reported in the summary output for the lm()<\/code> model.<\/p>\n
anova(mr) <\/code><\/pre>\n
##                 Analysis of Variance          Response: Species \r\n## \r\n##  Factor     d.f. Partial SS   MS           F     P     \r\n##  Area        1   4.237718e+03 4.237718e+03  1.14 0.2963\r\n##  Elevation   1   1.317666e+05 1.317666e+05 35.44 <.0001\r\n##  Nearest     1   2.797576e-01 2.797576e-01  0.00 0.9932\r\n##  Scruz       1   4.635787e+03 4.635787e+03  1.25 0.2752\r\n##  Adjacent    1   6.640639e+04 6.640639e+04 17.86 0.0003\r\n##  REGRESSION  5   2.918500e+05 5.837000e+04 15.70 <.0001\r\n##  ERROR      24   8.923137e+04 3.717974e+03<\/code><\/pre>\n
To run these same Partial F tests for the lm()<\/code> object we can use the base R drop1()<\/code> function.<\/p>\n
drop1(m, test = "F")<\/code><\/pre>\n
## Single term deletions\r\n## \r\n## Model:\r\n## Species ~ Area + Elevation + Nearest + Scruz + Adjacent\r\n##           Df Sum of Sq    RSS    AIC F value    Pr(>F)    \r\n## <none>                  89231 251.93                      \r\n## Area       1      4238  93469 251.33  1.1398 0.2963180    \r\n## Elevation  1    131767 220998 277.14 35.4404 3.823e-06 ***\r\n## Nearest    1         0  89232 249.93  0.0001 0.9931506    \r\n## Scruz      1      4636  93867 251.45  1.2469 0.2752082    \r\n## Adjacent   1     66406 155638 266.62 17.8609 0.0002971 ***\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1<\/code><\/pre>\n
The {rms} package offers a convenient plot method for anova()<\/code> objects that \u201cdraws dot charts depicting the importance of variables in the model\u201d (from the anova.rms<\/code> help page). The default importance measure is the chi-square statistic for each factor minus its degrees of freedom. There are several other importance measures available. See the anova.rms<\/code> help page for the what<\/code> argument.<\/p>\n
plot(anova(mr))<\/code><\/pre>\n
<\/p>\n<\/div>\n
\n
IQR effect plots<\/h2>\n
Earlier we noted that we don\u2019t use the summary()<\/code> function on ols()<\/code> model objects to see a model summary. That doesn\u2019t mean there isn\u2019t a summary()<\/code> method available. There is, but it does something entirely different than the summary method for lm()<\/code> objects. If we try it right now we get an error message:<\/p>\n
summary(mr)<\/code><\/pre>\n
## Error in summary.rms(mr): adjustment values not defined here or with datadist for Area Elevation Nearest Scruz Adjacent<\/code><\/pre>\n
Notice what it says: \u201cadjustment values not defined here or with datadist\u201d. This tells us we need to define \u201cadjustment values\u201d to use the summary()<\/code> function. Instead of summarizing the model, the {rms} summary function when called on {rms} model objects returns a \u201csummary of the effects of each factor\u201d. Let\u2019s demonstrate what this means.<\/p>\n
The easiest way to define adjustment values is to use the datadist()<\/code> function on our data frame and then set the global datadist option using the base R options()<\/code> function. Harrell frequently assigns datadist results to \u201cd\u201d so we do the same. (Also, the help page for datadist states, \u201cThe best method is probably to run datadist once before any models are fitted, storing the distribution summaries for all potential variables.\u201d I elected to wait for presentation purposes.)<\/p>\n
d <- datadist(gala)\r\noptions(datadist = "d")<\/code><\/pre>\n
If we print \u201cd\u201d we see adjustment levels for all variables in the model.<\/p>\n
d<\/code><\/pre>\n
##                 Species Endemics      Area Elevation Nearest   Scruz  Adjacent\r\n## Low:effect         13.0     7.25    0.2575     97.75   0.800  11.025    0.5200\r\n## Adjust to          42.0    18.00    2.5900    192.00   3.050  46.650    2.5900\r\n## High:effect        96.0    32.25   59.2375    435.25  10.025  81.075   59.2375\r\n## Low:prediction      2.0     1.45    0.0390     47.35   0.445   0.490    0.1000\r\n## High:prediction   319.1    85.40  782.6215   1229.40  40.205 193.905  782.6215\r\n## Low                 2.0     0.00    0.0100     25.00   0.200   0.000    0.0300\r\n## High              444.0    95.00 4669.3200   1707.00  47.400 290.200 4669.3200<\/code><\/pre>\n
Now let\u2019s call summary()<\/code> on our {rms} model object and see what it returns.<\/p>\n
summary(mr)<\/code><\/pre>\n
##              Effects              Response : Species \r\n## \r\n##  Factor    Low     High    Diff.   Effect     S.E.    Lower 0.95 Upper 0.95\r\n##  Area       0.2575  59.238  58.980  -1.411900  1.3225  -4.1413     1.3176  \r\n##  Elevation 97.7500 435.250 337.500 107.820000 18.1110  70.4400   145.2000  \r\n##  Nearest    0.8000  10.025   9.225   0.084353  9.7244 -19.9860    20.1550  \r\n##  Scruz     11.0250  81.075  70.050 -16.849000 15.0890 -47.9910    14.2930  \r\n##  Adjacent   0.5200  59.238  58.718  -4.392400  1.0393  -6.5374    -2.2473<\/code><\/pre>\n
This produces an interquartile range \u201cEffects\u201d summary for all predictors. For example, in the first row we see the effect of Area is about -1.41. To calculate this we predict Species when Area = 59.238 (High, 75th percentile) and when Area = 0.2575 (Low, 25th percentile) and take the difference in predicted species. For each prediction, all other variables are held at their \u201cAdjust to\u201d level as shown above when we printed the datadist object. In addition to the effect, the standard error (SE) of the effect and a 95% confidence interval on the effect is returned. In this case we\u2019re not sure if the effect is positive or negative. <\/p>\n
To calculate this effect measure using our lm()<\/code> model object we can use the predict()<\/code> function to make two predictions and then take the difference using the diff()<\/code> function. Notice all the values in the newdata argument come from the datadist object above. The two values for Area are the \u201cLow:effect\u201d and \u201cHigh:effect\u201d values. The other values from the \u201cAdjust to\u201d row.<\/p>\n
# How to get Area effect in summary(mr) output = -1.411900\r\np <- predict(m, newdata = data.frame(Area = c(0.2575, 59.238), \r\n                                Elevation=97.75, \r\n                                Nearest=3.05,\r\n                                Scruz=46.65,\r\n                                Adjacent=2.59))\r\np<\/code><\/pre>\n
##        1        2 \r\n## 26.90343 25.49153<\/code><\/pre>\n
diff(p)<\/code><\/pre>\n
##         2 \r\n## -1.411895<\/code><\/pre>\n
To get the standard error and confidence interval we can do the following. (Recall 58.980 is the IQR of Area.)<\/p>\n
se <- sqrt((58.980 * vcov(m)["Area","Area"] * 58.980))\r\nse<\/code><\/pre>\n
## [1] 1.32247<\/code><\/pre>\n
diff(p) + c(-1,1) * se * qt(0.975, df = 24)<\/code><\/pre>\n
## [1] -4.141340  1.317549<\/code><\/pre>\n
This is made a little easier using the emmeans()<\/code> function in the {emmeans} package. We use the at<\/code> argument to specify at what values we wish to make predictions for Area. All other values are held at their median by setting cov.reduce = median<\/code>. Sending the result to the {emmeans} function contrast<\/code> with argument \u201crevpairwise\u201d says to subtract the estimated means in reverse order. Finally we pipe into the confint()<\/code> to replicate the IQR effect estimate for Area that we saw in the ols()<\/code> summary.<\/p>\n
library(emmeans)\r\nemmeans(m, specs = "Area", at = list(Area = c(0.2575, 59.238)), \r\n        cov.reduce = median) |>\r\n  contrast("revpairwise") |>\r\n  confint()<\/code><\/pre>\n
##  contrast                estimate   SE df lower.CL upper.CL\r\n##  Area59.238 - Area0.2575    -1.41 1.32 24    -4.14     1.32\r\n## \r\n## Confidence level used: 0.95<\/code><\/pre>\n
There is also a plot method for {rms} summary objects. It plots the Effects and the 90, 95, and 99 percent confidence intervals using different shades of blue. Below we see the Area effect of -1.41 with relatively tight confidence intervals hovering around 0.<\/p>\n
plot(summary(mr))<\/code><\/pre>\n
<\/p>\n<\/div>\n
\n
Non-linear effects<\/h2>\n
Harrell advocates using non-linear effects in the form of regression splines. This makes a lot of sense when you pause and consider how many effects in real life are truly linear. Not many. Very few associations in nature indefinitely follow a straight line relationship. Fortunately, we can easily implement regression splines in R and specify how much non-linearity we want to entertain. We do this either in the form of degrees of freedom or knots. More of both means more non-linearity. I personally think of it as the number of times the relationship might change direction.<\/p>\n
One way to implement regression splines in R is via the ns()<\/code> function in the {splines} package, which comes installed with R. Below we fit a model that allows the effect of Nearest to change directions three times by specifying df=3<\/code>. We might think of this as sort of like using a 3-degree polynomial to model Nearest (but that\u2019s not what we\u2019re doing). The summary shows three coefficients for Nearest, neither of which have any interpretation.<\/p>\n
library(splines)\r\nm2 <- lm(Species ~ Area + Elevation + ns(Nearest, df = 3) + \r\n           Scruz + Adjacent, \r\n         data = gala)\r\nsummary(m2)<\/code><\/pre>\n
## \r\n## Call:\r\n## lm(formula = Species ~ Area + Elevation + ns(Nearest, df = 3) + \r\n##     Scruz + Adjacent, data = gala)\r\n## \r\n## Residuals:\r\n##     Min      1Q  Median      3Q     Max \r\n## -79.857 -28.285  -0.775  25.498 163.069 \r\n## \r\n## Coefficients:\r\n##                        Estimate Std. Error t value Pr(>|t|)    \r\n## (Intercept)            13.18825   27.39386   0.481   0.6350    \r\n## Area                   -0.03871    0.02098  -1.845   0.0785 .  \r\n## Elevation               0.35011    0.05086   6.883 6.51e-07 ***\r\n## ns(Nearest, df = 3)1 -168.44016   68.41261  -2.462   0.0221 *  \r\n## ns(Nearest, df = 3)2  -56.97623   57.67426  -0.988   0.3339    \r\n## ns(Nearest, df = 3)3   38.61120   46.22137   0.835   0.4125    \r\n## Scruz                  -0.13735    0.19900  -0.690   0.4973    \r\n## Adjacent               -0.08254    0.01743  -4.734   0.0001 ***\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\r\n## \r\n## Residual standard error: 55.11 on 22 degrees of freedom\r\n## Multiple R-squared:  0.8247, Adjusted R-squared:  0.7689 \r\n## F-statistic: 14.78 on 7 and 22 DF,  p-value: 5.403e-07<\/code><\/pre>\n
To decide whether we should keep this non-linear effect, we could use the anova()<\/code> function to compare this updated more complex model to the original model with only linear effects. The null of the test is that both models are equally adequate. It appears there is some evidence that this non-linearity improves the model.<\/p>\n
anova(m, m2)<\/code><\/pre>\n
## Analysis of Variance Table\r\n## \r\n## Model 1: Species ~ Area + Elevation + Nearest + Scruz + Adjacent\r\n## Model 2: Species ~ Area + Elevation + ns(Nearest, df = 3) + Scruz + Adjacent\r\n##   Res.Df   RSS Df Sum of Sq      F  Pr(>F)  \r\n## 1     24 89231                              \r\n## 2     22 66808  2     22424 3.6921 0.04144 *\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1<\/code><\/pre>\n
To entertain non-linear effects in {rms} use the rcs()<\/code> function, which stands for restricted cubic splines. Instead of degrees of freedom we specify knots. Three degrees of freedom corresponds to four knots.<\/p>\n
mr2 <- ols(Species ~ Area + Elevation + rcs(Nearest, 4) + \r\n               Scruz + Adjacent, \r\n             data = gala, x = TRUE)\r\nmr2<\/code><\/pre>\n
## Linear Regression Model\r\n## \r\n## ols(formula = Species ~ Area + Elevation + rcs(Nearest, 4) + \r\n##     Scruz + Adjacent, data = gala, x = TRUE)\r\n## \r\n##                  Model Likelihood    Discrimination    \r\n##                        Ratio Test           Indexes    \r\n## Obs       30    LR chi2     51.32    R2       0.819    \r\n## sigma55.9556    d.f.            7    R2 adj   0.762    \r\n## d.f.      22    Pr(> chi2) 0.0000    g      107.069    \r\n## \r\n## Residuals\r\n## \r\n##     Min      1Q  Median      3Q     Max \r\n## -90.884 -29.870  -3.715  28.019 163.497 \r\n## \r\n## \r\n##           Coef      S.E.      t     Pr(>|t|)\r\n## Intercept   14.1147   30.2230  0.47 0.6451  \r\n## Area        -0.0368    0.0212 -1.74 0.0964  \r\n## Elevation    0.3444    0.0512  6.73 <0.0001 \r\n## Nearest      3.6081   16.9795  0.21 0.8337  \r\n## Nearest'  -714.5143  891.8922 -0.80 0.4316  \r\n## Nearest'' 1001.4929 1207.4895  0.83 0.4158  \r\n## Scruz       -0.2286    0.1978 -1.16 0.2602  \r\n## Adjacent    -0.0806    0.0175 -4.60 0.0001<\/code><\/pre>\n
Notice we get three coefficients for Nearest that differ in value from our lm()<\/code> result. A brief explanation of why that is can be found on datamethods.org<\/a>. This paper<\/a> provides a deeper explanation.<\/p>\n
Calling the {rms} anova()<\/code> method on the model object returns a test for nonlinearity under the Nearest test, labeled \u201cNonlinear\u201d. Notice the resulting p-value is slightly higher and exceeds 0.05.<\/p>\n
anova(mr2)<\/code><\/pre>\n
##                 Analysis of Variance          Response: Species \r\n## \r\n##  Factor     d.f. Partial SS MS         F     P     \r\n##  Area        1     9446.295   9446.295  3.02 0.0964\r\n##  Elevation   1   141872.398 141872.398 45.31 <.0001\r\n##  Nearest     3    20348.931   6782.977  2.17 0.1208\r\n##   Nonlinear  2    20348.651  10174.326  3.25 0.0580\r\n##  Scruz       1     4182.177   4182.177  1.34 0.2602\r\n##  Adjacent    1    66204.433  66204.433 21.14 0.0001\r\n##  REGRESSION  7   312198.651  44599.807 14.24 <.0001\r\n##  ERROR      22    68882.715   3131.033<\/code><\/pre>\n
To replicate the {rms} anova()<\/code> output using lm()<\/code>, we need to change ns()<\/code> arguments. The following R code comes from this hbiostat.org page<\/a><\/p>\n
w <- rcs(gala$Nearest, 4)\r\nkn <- attr(w, 'parms')\r\nm2 <- lm(Species ~ Area + Elevation + ns(Nearest, knots = kn[2:3], \r\n                                         Boundary.knots = c(kn[1], kn[4])) \r\n         + Scruz + Adjacent, data = gala)\r\nanova(m, m2)<\/code><\/pre>\n
## Analysis of Variance Table\r\n## \r\n## Model 1: Species ~ Area + Elevation + Nearest + Scruz + Adjacent\r\n## Model 2: Species ~ Area + Elevation + ns(Nearest, knots = kn[2:3], Boundary.knots = c(kn[1], \r\n##     kn[4])) + Scruz + Adjacent\r\n##   Res.Df   RSS Df Sum of Sq      F  Pr(>F)  \r\n## 1     24 89231                              \r\n## 2     22 68883  2     20349 3.2495 0.05801 .\r\n## ---\r\n## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1<\/code><\/pre>\n
Which approach is better: ns()<\/code> or rcs()<\/code>? I don\u2019t pretend to know and I\u2019m not sure it really matters. The substance of the results doesn\u2019t differ much. Yes, we saw a p-value jump from 0.04 to 0.05, but we know better<\/a> than to make hard binary decisions based on a p-value.<\/p>\n<\/div>\n
\n
Effect plots<\/h2>\n
Since interpreting non-linear effect coefficients is all but impossible, we turn to visualization. When fitting a non-linear model with lm()<\/code>, the {effects} and {ggeffects} packages are fantastic for creating effect plots.<\/p>\n
m2 <- lm(Species ~ Area + Elevation + ns(Nearest, df = 3) + \r\n           Scruz + Adjacent, gala)\r\nlibrary(effects)\r\nEffect("Nearest", m2) |> plot()<\/code><\/pre>\n
<\/p>\n
library(ggeffects)\r\nggeffect(m2, "Nearest") |> plot()<\/code><\/pre>\n
<\/p>\n
In both cases we see the effect of Nearest is negative for values 0 – 20, but then seems to increase for values above 20.<\/p>\n
To get a similar plot for an {rms} model, we can use the Predict()<\/code> function (note the capital \u201cP\u201d; this is an {rms} function) and then pipe into plot()<\/code>. You can also pipe into ggplot()<\/code> and plotp()<\/code> to get ggplot2 and plotly plots, respectively.<\/p>\n
Predict(mr2, Nearest) |> plot() <\/code><\/pre>\n
<\/p>\n
To get all effect plots in one go, we can use the {effects} function allEffects()<\/code> with the lm()<\/code> model.<\/p>\n
allEffects(m2) |> plot()<\/code><\/pre>\n
<\/p>\n
To do the same with {rms}, we use the Predict()<\/code> function with no predictors specified. Notice the y-axis has the same limits for all plots.<\/p>\n
Predict(mr2) |> plot()<\/code><\/pre>\n
<\/p>\n<\/div>\n
\n
Diagnostics<\/h2>\n
A linear model makes several assumptions, two of which are constant variance of residuals and normality of residuals. The plot()<\/code> method for lm()<\/code> objects makes this easy to assess graphically.<\/p>\n
op <- par(mfrow = c(2,2))\r\nplot(m2)<\/code><\/pre>\n
<\/p>\n
par(op)<\/code><\/pre>\n
The left two plots assess constant variance using two different types of residuals. The upper right plot assess normality. The bottom right helps identify influential observations. The islands of Isabela and Santa Cruz stand out as either not being well fit by the model or unduly influencing the fit of the model.<\/p>\n
We can use the plot.lm()<\/code> method on the ols()<\/code> model, but only the first one (i.e., which = 1<\/code>).<\/p>\n
plot(mr2, which = 1)<\/code><\/pre>\n
<\/p>\n
To create a QQ plot to assess normality of residuals we need to extract the residuals using the residuals() function and plot with qqnorm()<\/code> and qqline()<\/code><\/p>\n
r <- residuals(mr2, type="student")\r\nqqnorm(r)\r\nqqline(r)<\/code><\/pre>\n
<\/p>\n
To check for possible influential observations for an lm()<\/code> model we can use the base R function influence.measures()<\/code> and its associated summary()<\/code> method. See ?influence.measures<\/code> for a breakdown of what is returned. The columns that begin \u201cdfb\u201d are the DFBETA statistics. DFBETAS indicate the effect that deleting each observation has on the estimates of the regression coefficients. That\u2019s why there\u2019s a DFBETA for every coefficient in the model.<\/p>\n
summary(influence.measures(m2))<\/code><\/pre>\n
## Potentially influential observations of\r\n##   lm(formula = Species ~ Area + Elevation + ns(Nearest, df = 3) +      Scruz + Adjacent, data = gala) :\r\n## \r\n##              dfb.1_ dfb.Area dfb.Elvt dfb.n(N,d=3)1 dfb.n(N,d=3)2 dfb.n(N,d=3)3\r\n## Darwin        -0.02   0.07    -0.08     0.05         -0.10         -0.10       \r\n## Fernandina     0.14   0.19    -0.15     0.02          0.00          0.03       \r\n## Genovesa      -0.07  -0.08     0.12     0.13         -0.15         -0.43       \r\n## Isabela        0.05 -18.76_*   4.55_*  -1.69_*       -1.12_*        0.94       \r\n## SanCristobal  -0.09  -0.11     0.18    -0.18          0.22          0.49       \r\n## SantaCruz      0.43  -1.03_*   1.58_*  -0.57         -0.96         -0.38       \r\n## Wolf           0.00   0.00     0.00     0.01          0.00          0.00       \r\n##              dfb.Scrz dfb.Adjc dffit    cov.r    cook.d   hat     \r\n## Darwin         0.48     0.00     0.62     2.36_*   0.05     0.48  \r\n## Fernandina    -0.07    -0.91    -1.51    27.71_*   0.30     0.95_*\r\n## Genovesa       0.14    -0.10    -0.47     3.16_*   0.03     0.57  \r\n## Isabela       -0.54    -0.78   -26.78_*   0.45    50.94_*   0.98_*\r\n## SanCristobal  -0.18    -0.12     0.58     2.56_*   0.04     0.50  \r\n## SantaCruz     -0.17    -1.16_*   2.33_*   0.01     0.36     0.21  \r\n## Wolf           0.03     0.00     0.04     2.21_*   0.00     0.34<\/code><\/pre>\n
The {rms} package offers the which.influence()<\/code> and show.influence()<\/code> tandem to identify observations (via DFBETAS) that effect regression coefficients with their removal. Below we set cutoff = 0.4<\/code>. An asterisk is placed next to a variable when any of the coefficients associated with that variable change by more than 0.4 standard errors upon removal of the observation. Below we see that removing the Isabela observation changes four of the coefficients. (Also, this is why we set x = TRUE<\/code> in our call to ols()<\/code>, so we could use these functions.) The values displayed are the observed values for each observation.<\/p>\n
w <- which.influence(mr2, cutoff = 0.4)\r\nshow.influence(w, gala)<\/code><\/pre>\n
##             Count     Area Elevation Nearest  Scruz Adjacent\r\n## Darwin          1     2.33       168    34.1 *290.2     2.85\r\n## Fernandina      1   634.49      1494     4.3   95.3 *4669.32\r\n## Isabela         4 *4669.32     *1707     0.7 * 28.1 * 634.49\r\n## Pinta           3 *  59.56     * 777    29.1  119.6 * 129.49\r\n## SanSalvador     3 * 572.33     * 906     0.2   19.8 *   4.89\r\n## SantaCruz       3 * 903.82     * 864     0.6    0.0 *   0.52\r\n## SantaFe         1    24.08       259   *16.5   16.5     0.52\r\n## SantaMaria      5 * 170.92     * 640   * 2.6   49.2 *   0.10<\/code><\/pre>\n<\/div>\n
\n
More information<\/h2>\n
Got a question? Harrell runs a discussion board<\/a> for the {rms} package.<\/p>\n
Chapter 14<\/a> of these course notes<\/a> by Thomas Love are also helpful.<\/p>\n
And here\u2019s a nice blog post called An Introduction to the Harrell\u201cverse\u201d: Predictive Modeling using the Hmisc and rms Packages<\/a> by Nicholas Ollberding.<\/p>\n
Eventually you may want to buy his book, Regression Modeling Strategies<\/a>. It\u2019s not cheap, but it\u2019s not outrageous either. Considering how much information, advice, code, examples, and references it contains, I think it\u2019s a bargain. It\u2019s a text that will provide many months, if not years, of study.<\/p>\n","protected":false},"excerpt":{"rendered":"
For a few years now I\u2019ve been chipping away at Regression Modeling Strategies by Frank Harrell. Between the exposition, the… Read more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[12],"tags":[],"class_list":["post-833","post","type-post","status-publish","format-standard","hentry","category-regression"],"_links":{"self":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/833","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/comments?post=833"}],"version-history":[{"count":9,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/833\/revisions"}],"predecessor-version":[{"id":877,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/833\/revisions\/877"}],"wp:attachment":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/media?parent=833"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/categories?post=833"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/tags?post=833"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}