Creation of classes and methods allows use of familiar generic functions such as plot() and summary() to do specific tasks for an object of a specific class. Here a new class “fitAndPlot” is defined. This class object contains a model fit part and a ggplot2 layer part. A print method and summary method are defined for the model fit part, and a plot method is defined for the ggplot2 layer part.
Chapter 9: Object-Oriented Programming, The Art of R Programming: http://www.nostarch.com/artofr.htm
Fractional polynomials etc: http://rpubs.com/kaz_yos/PolySpline3
## Load MASS for dataset
library(MASS)
## Load ggplot2 for plotting
library(ggplot2)
## Load gridExtra for fine tuning
library(gridExtra)
## Load splines for bs()
library(splines)
## Generalized Additive Models
library(mgcv)
## Load
data(mcycle)
## Plot raw data as a scatter plot (Used as a base plot)
plotBase <- ggplot(data = mcycle,
mapping = aes(x = times, y = accel)) +
labs(title = "Motor Cycle dataset") +
layer(geom = "point") +
scale_x_continuous(breaks = seq(from = 0, to = 100, by = 20)) +
theme_bw() +
theme(legend.key = element_blank())
plotBase
## Define a helper function for ggplot2
stat_fun <- function(fun, args, mapping) {
stat_function(fun = fun, args = args, mapping = mapping, n = 1001, size = 1.5, alpha = 1)
}
## Define prediction function
PredFun <- function(x, model) {
## Create data for prediction
newData <- data.frame(times = x)
## Predict for it
predict(object = model, newdata = newData)
}
## Define a "fit and plot" function (dependent on thresholds)
FitLayer <- function(formula, data = mcycle, fitFun = lm, ...) {
## Fit model
modelFit <- fitFun(formula = formula,
data = data, ...)
## Return a layer
layerFit <- stat_fun(fun = PredFun, args = list(model = modelFit), mapping = aes(x = times))
## Create a result object
out <- list(modelFit = modelFit, layerFit = layerFit)
## Define a class
class(out) <- "fitAndPlot"
## Return result
return(out)
}
## Define a plot method that creates a ggplot2 plot
plot.fitAndPlot <- function(obj, plotBase = plotBase) {
print(plotBase + obj$layerFit)
}
## Define a print method that returns the regression result
print.fitAndPlot <- function(obj) {
print(obj$modelFit)
}
## Define a summary method that returns the regression summary
summary.fitAndPlot <- function(obj) {
summary(obj$modelFit)
}
## Check methods available for the fitAndPlot class
methods(class = "fitAndPlot")
## [1] plot.fitAndPlot print.fitAndPlot summary.fitAndPlot
Linear spline example
## Create a fit object
fitObjLinearSplines <- FitLayer(accel ~ bs(times, df = 10, knots = NULL, degree = 1))
## Check object class
class(fitObjLinearSplines)
## [1] "fitAndPlot"
## Check list elements
names(fitObjLinearSplines)
## [1] "modelFit" "layerFit"
## Implicit print method by using the object name
fitObjLinearSplines
##
## Call:
## fitFun(formula = formula, data = data)
##
## Coefficients:
## (Intercept) bs(times, df = 10, knots = NULL, degree = 1)1
## -1.522 -0.452
## bs(times, df = 10, knots = NULL, degree = 1)2 bs(times, df = 10, knots = NULL, degree = 1)3
## -9.033 -48.703
## bs(times, df = 10, knots = NULL, degree = 1)4 bs(times, df = 10, knots = NULL, degree = 1)5
## -102.407 -115.993
## bs(times, df = 10, knots = NULL, degree = 1)6 bs(times, df = 10, knots = NULL, degree = 1)7
## -23.846 54.093
## bs(times, df = 10, knots = NULL, degree = 1)8 bs(times, df = 10, knots = NULL, degree = 1)9
## 10.693 1.502
## bs(times, df = 10, knots = NULL, degree = 1)10
## 0.968
## Summary method creates the summary of the model fit part of the object
summary(fitObjLinearSplines)
##
## Call:
## fitFun(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73.34 -11.89 -0.57 11.15 53.57
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.522 10.572 -0.14 0.88575
## bs(times, df = 10, knots = NULL, degree = 1)1 -0.452 15.166 -0.03 0.97628
## bs(times, df = 10, knots = NULL, degree = 1)2 -9.033 12.588 -0.72 0.47440
## bs(times, df = 10, knots = NULL, degree = 1)3 -48.703 12.627 -3.86 0.00018 ***
## bs(times, df = 10, knots = NULL, degree = 1)4 -102.407 13.344 -7.67 4.6e-12 ***
## bs(times, df = 10, knots = NULL, degree = 1)5 -115.993 14.010 -8.28 1.8e-13 ***
## bs(times, df = 10, knots = NULL, degree = 1)6 -23.846 12.873 -1.85 0.06638 .
## bs(times, df = 10, knots = NULL, degree = 1)7 54.093 14.267 3.79 0.00023 ***
## bs(times, df = 10, knots = NULL, degree = 1)8 10.693 13.164 0.81 0.41820
## bs(times, df = 10, knots = NULL, degree = 1)9 1.502 12.987 0.12 0.90809
## bs(times, df = 10, knots = NULL, degree = 1)10 0.968 15.413 0.06 0.95000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.4 on 122 degrees of freedom
## Multiple R-squared: 0.801, Adjusted R-squared: 0.785
## F-statistic: 49.2 on 10 and 122 DF, p-value: <2e-16
## Plot method takes the object and the plotBase to create
plot(fitObjLinearSplines, plotBase)
Generalized additive model example
## Create another model
fitObjGam <- FitLayer(accel ~ s(times, bs = "cr"), fitFun = gam, print = T)
## Print method
fitObjGam
##
## Family: gaussian
## Link function: identity
##
## Formula:
## accel ~ s(times, bs = "cr")
##
## Estimated degrees of freedom:
## 8.39 total = 9.39
##
## GCV score: 544.5
## Summary method
summary(fitObjGam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## accel ~ s(times, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -25.55 1.95 -13.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(times) 8.39 8.87 53.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.783 Deviance explained = 79.7%
## GCV score = 544.48 Scale est. = 506.04 n = 133
## Plot method
plot(fitObjGam, plotBase)
