library(tidyverse)
library(here)
library(sloop)
set.seed(36528)

Motivation

Some code from last week:

# some fake data
true_pi <- 0.8
true_mu <- c(0,   2)
true_sigma <- c(1,   0.5)
n <- 250
pop <- rbinom(n, size = 1, prob = 1 - true_pi) + 1
x <- rnorm(n, mean = true_mu[pop], sd = true_sigma[pop])

dmix <- function(x, theta){
  pi <- theta[1]
  mu_1 <- theta[2]
  mu_2 <- theta[3]
  sigma_1 <- theta[4]
  sigma_2 <- theta[5]
  pi * dnorm(x, mean = mu_1, sd = sigma_1) +
   (1 - pi) * dnorm(x, mean = mu_2, sd = sigma_2) 
}
nllhood_mix <- function(theta, x){
  d <- dmix(x, theta)
  -1 * sum(log(d))
}

Recall:

mle1 <- optim(par = c(1, 0, 1.8, 1.2, 0.7), fn = nllhood_mix, x = x, method = "BFGS")
mle1
## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
## 
## $value
## [1] 411.6919
## 
## $counts
## function gradient 
##      125       41 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
mle2 <- optim(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")
mle2
## $par
## [1] -3410.0610 -6066.8409  -230.4462   550.6597   588.7451
## 
## $value
## [1] -190.2385
## 
## $counts
## function gradient 
##      185      100 
## 
## $convergence
## [1] 1
## 
## $message
## NULL

Your Turn What do you do with this output? How might you improve what the output looks like to help you as a user?

Our goal: To get a more user friendly output from optim()

Object Oriented Programming (OOP)

A programming paradigm where objects are central. Two important terms:

Gives us a way to have the same function operate differently on different classes of object, e.g. summary()

z <- factor(sample(letters[1:5], size = 10, replace = TRUE))
y <- rnorm(n = 10)
summary(z)
## a b c d e 
## 1 2 3 3 1
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.6029 -0.3881  0.3527  0.2827  1.0383  2.3370

There is more than one OOP system in R, we’ll focus on S3.

S3

The most common OO system R. Also the simplest.

S3 allows your functions to return rich results with user-friendly display and programmer-friendly internals. S3 is used throughout base R, so it’s important to master if you want to extend base R functions to work with new types of input.

Advanced R

Why bother?

Primarily to make your functions easier to use (for you or others), by creating a new class and defining some methods:

  • Have your functions print nice output
  • Implement methods for familiar functions, e.g. plot(), summary()

More advanced make it easy for others to extend your work for other data types, by supplying generic functions.

Orientation

As our functions get more complicated they need to return complicated things.

E.g. Linear models

mod <- lm(mpg ~ wt, data = mtcars)
str(mod)
## List of 12
##  $ coefficients : Named num [1:2] 37.29 -5.34
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "wt"
##  $ residuals    : Named num [1:32] -2.28 -0.92 -2.09 1.3 -0.2 ...
##   ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
##  $ effects      : Named num [1:32] -113.65 -29.116 -1.661 1.631 0.111 ...
##   ..- attr(*, "names")= chr [1:32] "(Intercept)" "wt" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:32] 23.3 21.9 24.9 20.1 18.9 ...
##   ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:32, 1:2] -5.657 0.177 0.177 0.177 0.177 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "wt"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.18 1.05
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 30
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = mpg ~ wt, data = mtcars)
##  $ terms        :Classes 'terms', 'formula'  language mpg ~ wt
##   .. ..- attr(*, "variables")= language list(mpg, wt)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "mpg" "wt"
##   .. .. .. ..$ : chr "wt"
##   .. ..- attr(*, "term.labels")= chr "wt"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(mpg, wt)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "mpg" "wt"
##  $ model        :'data.frame':   32 obs. of  2 variables:
##   ..$ mpg: num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##   ..$ wt : num [1:32] 2.62 2.88 2.32 3.21 3.44 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language mpg ~ wt
##   .. .. ..- attr(*, "variables")= language list(mpg, wt)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "mpg" "wt"
##   .. .. .. .. ..$ : chr "wt"
##   .. .. ..- attr(*, "term.labels")= chr "wt"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(mpg, wt)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "mpg" "wt"
##  - attr(*, "class")= chr "lm"

Complicated objects are built on top of simple objects (vectors, matrices, lists etc.) with the additional of attributes.

E.g. lm() returns a list along with the a class attribute.

The class attribute is special, it’s how R knows what kind of object something is.

Use class() to see the class,

class(mod) # the class of this object
## [1] "lm"

typeof() to see the base type the object is built on,

typeof(mod) # the base type this object is built on
## [1] "list"

and attributes() to see all attributes:

attributes(mod)
## $names
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"        
## 
## $class
## [1] "lm"

Your turn What class is z? What base object is it built on? What other attributes does this object have?

z <- factor(sample(letters[1:5], size = 10, replace = TRUE))

class(z)
## [1] "factor"
typeof(z)
## [1] "integer"
attributes(z)
## $levels
## [1] "a" "b" "c" "d" "e"
## 
## $class
## [1] "factor"

Does everything have a class? class() will return something, even if the object doesn’t have a class attribute:

class(y)
## [1] "numeric"
attributes(y)
## NULL

Methods

In S3 methods are associated with generic functions (generics for short).

This is different to most other OO languages out there, which use encapsulated OO - the methods are associated with the object.

When an S3 generic function is called, it looks at the class of the object passed to it, and searches for the appropriate method, of the form generic.class, e.g.

summary(mod) # calls summary.lm
## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10
class(mod)
## [1] "lm"

looks for finds and calls summary.lm()

This process is known as method dispatch. You can see the process with s3_dispatch() in the sloop package:

s3_dispatch(summary(mod))
## => summary.lm
##  * summary.default

=> the method that is called, * a method that is defined but not called.

If you want to know if a function is a generic you can use ftype() in sloop

ftype(summary)
## [1] "S3"      "generic"
ftype(summary.lm)
## [1] "S3"     "method"

Your turn You can use . in a name in R without it being a method. Can you figure out if t.test() is the t method for objects of class test? What about t.data.frame()?

# install.packages("sloop")
library(sloop)

ftype(t.test) # S3 generic, not a method for `t()` for class "test"
## [1] "S3"      "generic"
ftype(t.data.frame) # S3 method for `t()` for class "data.frame"
## [1] "S3"     "method"
fake <- data.frame(x = 1:10)
class(fake)
## [1] "data.frame"
s3_dispatch(t(fake))  # let's you see how R found the method
## => t.data.frame
## -> t.default

Unanswered question: What if we define a frame class and a t.data() generic, then and a frame method for t.data(), do things break?

class doesn’t get any special treatment

You can set the class of an object just like any usual assignment:

class(mod) <- "table"
mod
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        coefficients 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                37.285126, -5.344472 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           residuals 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     -2.2826106, -0.9197704, -2.0859521, 1.2973499, -0.2001440, -0.6932545, -3.9053627, 4.1637381, 2.3499593, 0.2998560, -1.1001440, 0.8668731, -0.0502472, -1.8830236, 1.1733496, 2.1032876, 5.9810744, 6.8727113, 1.7461954, 6.4219792, -2.6110037, -2.9725862, -3.7268663, -3.4623553, 2.4643670, 0.3564263, 0.1520430, 1.2010593, -4.5431513, -2.7809399, -3.2053627, -1.0274952 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             effects 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    -113.6497374, -29.1157217, -1.6613339, 1.6313943, 0.1111305, -0.3840041, -3.6072442, 4.5003125, 2.6905817, 0.6111305, -0.7888695, 1.1143917, 0.2316793, -1.6061571, 1.3014525, 2.2137818, 6.0995633, 7.3094734, 2.2421594, 6.8956792, -2.2010595, -2.6694078, -3.4150859, -3.1915608, 2.7346556, 0.8200064, 0.5948771, 1.7073457, -4.2045529, -2.4018616, -2.9072442, -0.6494289 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                rank 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   2 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       fitted.values 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         23.282611, 21.919770, 24.885952, 20.102650, 18.900144, 18.793255, 18.205363, 20.236262, 20.450041, 18.900144, 18.900144, 15.533127, 17.350247, 17.083024, 9.226650, 8.296712, 8.718926, 25.527289, 28.653805, 27.478021, 24.111004, 18.472586, 18.926866, 16.762355, 16.735633, 26.943574, 25.847957, 29.198941, 20.343151, 22.480940, 18.205363, 22.427495 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              assign 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                0, 1 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  qr 
## -5.656854249, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, -18.199514334, 5.447820482, 0.148230003, -0.016055881, -0.057356801, -0.061027994, -0.081219555, -0.011466889, -0.004124504, -0.057356801, -0.057356801, -0.172999378, -0.110589098, -0.119767081, -0.389599760, -0.421539139, -0.407037927, 0.170257160, 0.277639553, 0.237256431, 0.121613854, -0.072041573, -0.056439003, -0.130780659, -0.131698458, 0.218900467, 0.181270739, 0.296362637, -0.007795696, 0.065628162, -0.081219555, 0.063792566, 1.176776695, 1.046354399, 1.000000000, 2.000000000, 0.000000100, 2.000000000 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         df.residual 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  30 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             xlevels 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                NULL 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                call 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               lm(formula = mpg ~ wt, data = mtcars) 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               terms 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            mpg ~ wt 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               model 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                      21.000, 21.000, 22.800, 21.400, 18.700, 18.100, 14.300, 24.400, 22.800, 19.200, 17.800, 16.400, 17.300, 15.200, 10.400, 10.400, 14.700, 32.400, 30.400, 33.900, 21.500, 15.500, 15.200, 13.300, 19.200, 27.300, 26.000, 30.400, 15.800, 19.700, 15.000, 21.400, 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3.150, 3.440, 3.440, 4.070, 3.730, 3.780, 5.250, 5.424, 5.345, 2.200, 1.615, 1.835, 2.465, 3.520, 3.435, 3.840, 3.845, 1.935, 2.140, 1.513, 3.170, 2.770, 3.570, 2.780

But you probably wouldn’t want to…

If we remove the class of mod

class(mod) <- NULL
mod 
## $coefficients
## (Intercept)          wt 
##   37.285126   -5.344472 
## 
## $residuals
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##          -2.2826106          -0.9197704          -2.0859521 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##           1.2973499          -0.2001440          -0.6932545 
##          Duster 360           Merc 240D            Merc 230 
##          -3.9053627           4.1637381           2.3499593 
##            Merc 280           Merc 280C          Merc 450SE 
##           0.2998560          -1.1001440           0.8668731 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##          -0.0502472          -1.8830236           1.1733496 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##           2.1032876           5.9810744           6.8727113 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##           1.7461954           6.4219792          -2.6110037 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##          -2.9725862          -3.7268663          -3.4623553 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##           2.4643670           0.3564263           0.1520430 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##           1.2010593          -4.5431513          -2.7809399 
##       Maserati Bora          Volvo 142E 
##          -3.2053627          -1.0274952 
## 
## $effects
##  (Intercept)           wt                                        
## -113.6497374  -29.1157217   -1.6613339    1.6313943    0.1111305 
##                                                                  
##   -0.3840041   -3.6072442    4.5003125    2.6905817    0.6111305 
##                                                                  
##   -0.7888695    1.1143917    0.2316793   -1.6061571    1.3014525 
##                                                                  
##    2.2137818    6.0995633    7.3094734    2.2421594    6.8956792 
##                                                                  
##   -2.2010595   -2.6694078   -3.4150859   -3.1915608    2.7346556 
##                                                                  
##    0.8200064    0.5948771    1.7073457   -4.2045529   -2.4018616 
##                           
##   -2.9072442   -0.6494289 
## 
## $rank
## [1] 2
## 
## $fitted.values
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##           23.282611           21.919770           24.885952 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##           20.102650           18.900144           18.793255 
##          Duster 360           Merc 240D            Merc 230 
##           18.205363           20.236262           20.450041 
##            Merc 280           Merc 280C          Merc 450SE 
##           18.900144           18.900144           15.533127 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##           17.350247           17.083024            9.226650 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##            8.296712            8.718926           25.527289 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##           28.653805           27.478021           24.111004 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##           18.472586           18.926866           16.762355 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##           16.735633           26.943574           25.847957 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##           29.198941           20.343151           22.480940 
##       Maserati Bora          Volvo 142E 
##           18.205363           22.427495 
## 
## $assign
## [1] 0 1
## 
## $qr
## $qr
##                     (Intercept)            wt
## Mazda RX4            -5.6568542 -18.199514334
## Mazda RX4 Wag         0.1767767   5.447820482
## Datsun 710            0.1767767   0.148230003
## Hornet 4 Drive        0.1767767  -0.016055881
## Hornet Sportabout     0.1767767  -0.057356801
## Valiant               0.1767767  -0.061027994
## Duster 360            0.1767767  -0.081219555
## Merc 240D             0.1767767  -0.011466889
## Merc 230              0.1767767  -0.004124504
## Merc 280              0.1767767  -0.057356801
## Merc 280C             0.1767767  -0.057356801
## Merc 450SE            0.1767767  -0.172999378
## Merc 450SL            0.1767767  -0.110589098
## Merc 450SLC           0.1767767  -0.119767081
## Cadillac Fleetwood    0.1767767  -0.389599760
## Lincoln Continental   0.1767767  -0.421539139
## Chrysler Imperial     0.1767767  -0.407037927
## Fiat 128              0.1767767   0.170257160
## Honda Civic           0.1767767   0.277639553
## Toyota Corolla        0.1767767   0.237256431
## Toyota Corona         0.1767767   0.121613854
## Dodge Challenger      0.1767767  -0.072041573
## AMC Javelin           0.1767767  -0.056439003
## Camaro Z28            0.1767767  -0.130780659
## Pontiac Firebird      0.1767767  -0.131698458
## Fiat X1-9             0.1767767   0.218900467
## Porsche 914-2         0.1767767   0.181270739
## Lotus Europa          0.1767767   0.296362637
## Ford Pantera L        0.1767767  -0.007795696
## Ferrari Dino          0.1767767   0.065628162
## Maserati Bora         0.1767767  -0.081219555
## Volvo 142E            0.1767767   0.063792566
## attr(,"assign")
## [1] 0 1
## 
## $qraux
## [1] 1.176777 1.046354
## 
## $pivot
## [1] 1 2
## 
## $tol
## [1] 1e-07
## 
## $rank
## [1] 2
## 
## attr(,"class")
## [1] "qr"
## 
## $df.residual
## [1] 30
## 
## $xlevels
## named list()
## 
## $call
## lm(formula = mpg ~ wt, data = mtcars)
## 
## $terms
## mpg ~ wt
## attr(,"variables")
## list(mpg, wt)
## attr(,"factors")
##     wt
## mpg  0
## wt   1
## attr(,"term.labels")
## [1] "wt"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(mpg, wt)
## attr(,"dataClasses")
##       mpg        wt 
## "numeric" "numeric" 
## 
## $model
##                      mpg    wt
## Mazda RX4           21.0 2.620
## Mazda RX4 Wag       21.0 2.875
## Datsun 710          22.8 2.320
## Hornet 4 Drive      21.4 3.215
## Hornet Sportabout   18.7 3.440
## Valiant             18.1 3.460
## Duster 360          14.3 3.570
## Merc 240D           24.4 3.190
## Merc 230            22.8 3.150
## Merc 280            19.2 3.440
## Merc 280C           17.8 3.440
## Merc 450SE          16.4 4.070
## Merc 450SL          17.3 3.730
## Merc 450SLC         15.2 3.780
## Cadillac Fleetwood  10.4 5.250
## Lincoln Continental 10.4 5.424
## Chrysler Imperial   14.7 5.345
## Fiat 128            32.4 2.200
## Honda Civic         30.4 1.615
## Toyota Corolla      33.9 1.835
## Toyota Corona       21.5 2.465
## Dodge Challenger    15.5 3.520
## AMC Javelin         15.2 3.435
## Camaro Z28          13.3 3.840
## Pontiac Firebird    19.2 3.845
## Fiat X1-9           27.3 1.935
## Porsche 914-2       26.0 2.140
## Lotus Europa        30.4 1.513
## Ford Pantera L      15.8 3.170
## Ferrari Dino        19.7 2.770
## Maserati Bora       15.0 3.570
## Volvo 142E          21.4 2.780
s3_dispatch(print(mod))
##    print.list
## => print.default

It just prints like any list. This isn’t very useful output…why?

To build a new object

You can either make the object from it’s base type and add a class (or any other attributes):

y <- 1:6
class(y) <- "my_special_kind_of_vector"
attr(y, "max") <- 6
y
## [1] 1 2 3 4 5 6
## attr(,"class")
## [1] "my_special_kind_of_vector"
## attr(,"max")
## [1] 6

Or use structure() to do it in one go

y <- structure(1:6, max = 6, class = "my_special_kind_of_vector")
y
## [1] 1 2 3 4 5 6
## attr(,"max")
## [1] 6
## attr(,"class")
## [1] "my_special_kind_of_vector"

Back to our goal

Let’s build a new class to hold optim() output, then write a print() method that prints output in more useful fashion.

General procedure:

  1. Decide on a name for the class
  2. Decide what structure/properties the class should have
  3. Write a low level constructor function creates this kind of object, and a user friendly helper function to easily get an example.
  4. Write methods

Your turn What shall we call the class to hold an optim() result?

  1. my_new_class, optim, my_optim, cl_optim, optim_class, optim_out, optimum

  2. Takes the same structure as existing optim() output - list with certain elements. We aren’t going to check for these elements in our constructor.

A constructor

https://adv-r.hadley.nz/s3.html#s3-constrcutor

First we’ll write a function that creates an object of class optimum.

By convention, call it new_classname():

new_optimum <- function(l){
  # This is where you need to enforce any desired structure
  stopifnot(is.list(l))
  
  structure(
    l,
    class = "optimum"
  )
}

new_optimum(list(par = c(1, 2, 3), y = "weird"))
## $par
## [1] 1 2 3
## 
## $y
## [1] "weird"
## 
## attr(,"class")
## [1] "optimum"

I’m being a little lazy checking the input, but I’ll assume the incoming list has all the components of an optim() output.

new_optimum(mle1)
## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
## 
## $value
## [1] 411.6919
## 
## $counts
## function gradient 
##      125       41 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## attr(,"class")
## [1] "optimum"
s3_dispatch(print(new_optimum(mle1)))
##    print.optimum
## => print.default

A helper

https://adv-r.hadley.nz/s3.html#helpers

It’s useful to have a function class_name() that conveniently creates object of this class, which we’ll assume just passes it’s arguments onto optim():

optimum <- function(..., quiet = TRUE){
  if(quiet){
    suppressWarnings({
      opt <- optim(...)
      })
  } else {
    opt <- optim(...)
  }
  new_optimum(opt)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7), 
  fn = nllhood_mix, x = x, method = "BFGS")  
## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
## 
## $value
## [1] 411.6919
## 
## $counts
## function gradient 
##      125       41 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## attr(,"class")
## [1] "optimum"

Our first method print.optimum()

To define a new method, you just write a regular function with a special name generic.class():

print.optimum <- function(x, ...){
  cat("Final parameter values: \n")
  cat(x$par)
  invisible(x)
}

optimum(par = c(1, 0, 1.8, 1.2, 0.7), 
  fn = nllhood_mix, x = x, method = "BFGS")  
## Final parameter values: 
## 0.9362468 0.3054438 2.248824 1.24127 0.1166007

Your turn Edit the method to just print the final parameter values.

Better output

We should check that value of convergence and print a message if optim() didn’t converge:

print.optimum <- function(x, ...){
  if(x$convergence != 0){
    cat("Failure to converge \n")
  } else {
    cat("Coverged! \n")
  }
  cat("Final parameter values \n")
  print(x$par)
  invisible(x)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7), 
  fn = nllhood_mix, x = x, method = "BFGS")  
## Coverged! 
## Final parameter values 
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
optimum(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")
## Failure to converge 
## Final parameter values 
## [1] -3410.0610 -6066.8409  -230.4462   550.6597   588.7451

Even better we might check some more convergence values and add some color! (might need install.packages("crayon")):

print.optimum <- function(x, ...){
  if(x$convergence != 0){
    cat(crayon::red("Failure to converge \n"))
    if (x$convergence == 1){
      cat("Iteration limit `maxit` reached. \n")
    } else {
      cat(x$message)
    }
  } else {
    cat(crayon::green("Coverged! \n"))
  }
  cat("Final parameter values \n")
  print(x$par)
  invisible(x)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7), 
  fn = nllhood_mix, x = x, method = "BFGS")  
## Coverged! 
## Final parameter values 
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
optimum(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")
## Failure to converge 
## Iteration limit `maxit` reached. 
## Final parameter values 
## [1] -3410.0610 -6066.8409  -230.4462   550.6597   588.7451

Another method

Your turn Is coef() a generic? Write a coef() method for optimum that returns the parameter estimates.

ftype(coef) # Yes, it's a generic
## [1] "S3"      "generic"
coef # Our method needs to accept the same arguments as the generic
## function (object, ...) 
## UseMethod("coef")
## <bytecode: 0x7fd1a7440f08>
## <environment: namespace:stats>
coef.optimum <- function(object, ...) {
  object$par
}

mle_1a <- optimum(par = c(1, 0, 1.8, 1.2, 0.7), 
  fn = nllhood_mix, x = x, method = "BFGS")  

coef(mle_1a)
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007

FYI stats4::mle()

Does a similar thing to what we have implemented but it uses the S4 OO system.