Skip to contents

There is a rich selection of R packages implementing algorithms for classification and regression tasks out there. The authors legitimately take the liberty to tailor the function interfaces according to their own taste and needs. For us other users, however, this often results in struggling with user interfaces, some of which are rather weird - to put it mildly - and almost always different in terms of arguments and result structures. ModTools pursues the goal of offering uniform handling for the most important regression and classification models in applied data analyses.
The function FitMod() is designed as a simple and consistent interface to these original functions while maintaining the flexibility to pass on all possible arguments. print, plot, summary and predict operations can so be carried out following the same logic. The results will again be reshaped to a reasonable standard.

For all the functions of this package Google styleguides are used as naming rules (in absence of convincing alternatives). The 'BigCamelCase' style has been consequently applied to functions borrowed from contributed R packages as well.

As always: Feedback, feature requests, bugreports and other suggestions are welcome!

Warning

This package is still under development. You should be aware that everything in the package might be subject to change. Backward compatibility is not yet guaranteed. Functions may be deleted or renamed and new syntax may be inconsistent with earlier versions. By release of version 1.0 the "deprecated-defunct process" will be installed.

Details

The ModTools::FitMod()) function comprises interfaces to the following models:

Regression:
lm()Linear model OLS (base)
lmrob()Robust linear model (robustbase)
poisson()GLM model with family poisson (base)
negbin()GLM model with family negative.binomial (MASS)
gamma()GLM model with family gamma (base)
tobit()Tobit model for censored responses (package AER)
Classification:
lda()
Linear discriminant analysis (MASS)qda()
Quadratic discriminant analysis (MASS)logit()
Logistic Regression model glm, family binomial(logit)(base)multinom()
Multinomial Regression model (nnet)polr()
Proportional odds model (MASS)rpart()
Regression and classification trees (rpart)nnet()
Neuronal networks (nnet)randomForest()
Random forests (randomForest)C5.0()
C5.0 tree (C50)svm()
Support vector machines (e1071)naive_bayes()
Naive Bayes classificator (naivebayes)LogitBoost()
Logit boost (using decision stumps as weak learners) (ModTools)
Preprocess:
SplitTrainTest()Splits a data frame or index vector into a training and a test sample
OverSample()Get balanced datasets by sampling with replacement.
Manipulating rpart objects:
CP()
Extract and plot complexity table of an rpart tree.Node()
Accessor to the most important properties of a node, being a split or a leaf.Rules()
Extract the decision rules from top to the end node of an rpart tree.LeafRates()
Returns the misclassification rates in all end nodes.
Prediction and Validation:
Response()Extract the response variable of any model.
predict()Consistent predict for FitMod models
VarImp()Variable importance for most FitMod models
ROC()ROC curves for all dichotomous classification FitMod models
BestCut()Find the optimal cut for a classification based on the ROC curve.
PlotLift()Produces a lift chart for a binary classification model
TModC()Aggregated results for multiple FitMod classification models
Tune()Tuning approaches to find optimal parameters for FitMod classification models.
RobSummary()Robust summary for GLM models (poisson).
Tests:
BreuschPaganTest()
Breusch-Pagan test against heteroskedasticity.

Author

Andri Signorell
Helsana Versicherungen AG, Health Sciences, Zurich
HWZ University of Applied Sciences in Business Administration Zurich.

Includes R source code and/or documentation previously published by (in alphabetical order):
Bernhard Compton, Marcel Dettling, Max Kuhn, Michal Majka, Dan Putler, Jarek Tuszynski, Robin Xavier, Achim Zeileis

The good things come from all these guys, any problems are likely due to my tweaking. Thank you all!

Maintainer: Andri Signorell <andri@signorell.net>

Examples


r.swiss <- FitMod(Fertility ~ ., swiss, fitfn="lm")
r.swiss
#> 
#> Call:
#> lm(formula = Fertility ~ ., data = swiss)
#> 
#> Coefficients:
#>      (Intercept)       Agriculture       Examination         Education  
#>          66.9152           -0.1721           -0.2580           -0.8709  
#>         Catholic  Infant.Mortality  
#>           0.1041            1.0770  
#> 
# PlotTA(r.swiss)
# PlotQQNorm(r.swiss)


## Count models

data(housing, package="MASS")

# poisson count
r.pois <- FitMod(Freq ~ Infl*Type*Cont + Sat, family=poisson, data=housing, fitfn="poisson")

# negative binomial count
r.nb <- FitMod(Freq ~ Infl*Type*Cont + Sat, data=housing, fitfn="negbin")
summary(r.nb)
#> 
#> Call:
#> glm.nb(formula = Freq ~ Infl * Type * Cont + Sat, data = housing, 
#>     init.theta = 11.58904011, link = log)
#> 
#> Coefficients:
#>                                   Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)                        3.14037    0.20801  15.097  < 2e-16 ***
#> InflMedium                         0.26615    0.28836   0.923 0.356028    
#> InflHigh                          -0.25336    0.30111  -0.841 0.400113    
#> TypeApartment                      0.39532    0.28603   1.382 0.166945    
#> TypeAtrium                        -0.77634    0.32149  -2.415 0.015742 *  
#> TypeTerrace                       -0.80653    0.32296  -2.497 0.012515 *  
#> ContHigh                          -0.02395    0.29474  -0.081 0.935244    
#> Sat.L                              0.16860    0.07525   2.240 0.025070 *  
#> Sat.Q                              0.25749    0.07821   3.292 0.000994 ***
#> InflMedium:TypeApartment          -0.12807    0.39875  -0.321 0.748082    
#> InflHigh:TypeApartment             0.15140    0.41087   0.368 0.712507    
#> InflMedium:TypeAtrium             -0.41265    0.45633  -0.904 0.365849    
#> InflHigh:TypeAtrium               -0.13006    0.47471  -0.274 0.784108    
#> InflMedium:TypeTerrace             0.01666    0.44423   0.038 0.970080    
#> InflHigh:TypeTerrace              -0.06969    0.47429  -0.147 0.883177    
#> InflMedium:ContHigh               -0.12632    0.41042  -0.308 0.758256    
#> InflHigh:ContHigh                 -0.59072    0.44424  -1.330 0.183609    
#> TypeApartment:ContHigh             0.52179    0.40019   1.304 0.192281    
#> TypeAtrium:ContHigh                0.71556    0.43760   1.635 0.102013    
#> TypeTerrace:ContHigh               1.15028    0.43270   2.658 0.007852 ** 
#> InflMedium:TypeApartment:ContHigh  0.01178    0.56001   0.021 0.983224    
#> InflHigh:TypeApartment:ContHigh    0.13235    0.59111   0.224 0.822837    
#> InflMedium:TypeAtrium:ContHigh     0.14135    0.62032   0.228 0.819752    
#> InflHigh:TypeAtrium:ContHigh       0.42364    0.65770   0.644 0.519496    
#> InflMedium:TypeTerrace:ContHigh   -0.51739    0.60483  -0.855 0.392307    
#> InflHigh:TypeTerrace:ContHigh     -0.49705    0.66522  -0.747 0.454946    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for Negative Binomial(11.589) family taken to be 1)
#> 
#>     Null deviance: 275.558  on 71  degrees of freedom
#> Residual deviance:  67.887  on 46  degrees of freedom
#> AIC: 535.29
#> 
#> Number of Fisher Scoring iterations: 1
#> 
#> 
#>               Theta:  11.59 
#>           Std. Err.:  2.92 
#> 
#>  2 x log-likelihood:  -481.289 

r.log <- FitMod(log(Freq) ~ Infl*Type*Cont + Sat, data=housing, fitfn="lm")
summary(r.log)
#> 
#> Call:
#> lm(formula = log(Freq) ~ Infl * Type * Cont + Sat, data = housing)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.90539 -0.22284  0.01119  0.18857  0.96514 
#> 
#> Coefficients:
#>                                    Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                        3.140416   0.267518  11.739 1.95e-15 ***
#> InflMedium                         0.259891   0.378328   0.687   0.4956    
#> InflHigh                          -0.379083   0.378328  -1.002   0.3216    
#> TypeApartment                      0.219444   0.378328   0.580   0.5647    
#> TypeAtrium                        -0.785497   0.378328  -2.076   0.0435 *  
#> TypeTerrace                       -0.931069   0.378328  -2.461   0.0177 *  
#> ContHigh                          -0.075612   0.378328  -0.200   0.8425    
#> Sat.L                              0.199559   0.094582   2.110   0.0403 *  
#> Sat.Q                              0.199376   0.094582   2.108   0.0405 *  
#> InflMedium:TypeApartment           0.048725   0.535036   0.091   0.9278    
#> InflHigh:TypeApartment             0.398374   0.535036   0.745   0.4603    
#> InflMedium:TypeAtrium             -0.400214   0.535036  -0.748   0.4583    
#> InflHigh:TypeAtrium                0.002462   0.535036   0.005   0.9963    
#> InflMedium:TypeTerrace             0.143412   0.535036   0.268   0.7899    
#> InflHigh:TypeTerrace               0.154151   0.535036   0.288   0.7746    
#> InflMedium:ContHigh               -0.105500   0.535036  -0.197   0.8446    
#> InflHigh:ContHigh                 -0.737873   0.535036  -1.379   0.1745    
#> TypeApartment:ContHigh             0.697935   0.535036   1.304   0.1986    
#> TypeAtrium:ContHigh                0.763012   0.535036   1.426   0.1606    
#> TypeTerrace:ContHigh               1.114096   0.535036   2.082   0.0429 *  
#> InflMedium:TypeApartment:ContHigh -0.141229   0.756655  -0.187   0.8528    
#> InflHigh:TypeApartment:ContHigh    0.087753   0.756655   0.116   0.9082    
#> InflMedium:TypeAtrium:ContHigh     0.060731   0.756655   0.080   0.9364    
#> InflHigh:TypeAtrium:ContHigh       0.503181   0.756655   0.665   0.5094    
#> InflMedium:TypeTerrace:ContHigh   -0.531148   0.756655  -0.702   0.4862    
#> InflHigh:TypeTerrace:ContHigh     -0.296310   0.756655  -0.392   0.6972    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 0.4634 on 46 degrees of freedom
#> Multiple R-squared:  0.751,	Adjusted R-squared:  0.6157 
#> F-statistic: 5.551 on 25 and 46 DF,  p-value: 2.796e-07
#> 

r.ols <- FitMod(Freq ~ Infl*Type*Cont + Sat, data=housing, fitfn="lm")
summary(r.ols)
#> 
#> Call:
#> lm(formula = Freq ~ Infl * Type * Cont + Sat, data = housing)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -22.4861  -5.2361  -0.1944   3.5347  27.0556 
#> 
#> Coefficients:
#>                                     Estimate Std. Error t value Pr(>|t|)   
#> (Intercept)                        2.333e+01  7.076e+00   3.298  0.00189 **
#> InflMedium                         7.333e+00  1.001e+01   0.733  0.46737   
#> InflHigh                          -4.333e+00  1.001e+01  -0.433  0.66700   
#> TypeApartment                      1.033e+01  1.001e+01   1.033  0.30717   
#> TypeAtrium                        -1.267e+01  1.001e+01  -1.266  0.21195   
#> TypeTerrace                       -1.300e+01  1.001e+01  -1.299  0.20037   
#> ContHigh                          -1.680e-14  1.001e+01   0.000  1.00000   
#> Sat.L                              2.976e+00  2.502e+00   1.190  0.24034   
#> Sat.Q                              5.835e+00  2.502e+00   2.332  0.02412 * 
#> InflMedium:TypeApartment          -1.667e+00  1.415e+01  -0.118  0.90676   
#> InflHigh:TypeApartment             3.333e+00  1.415e+01   0.236  0.81483   
#> InflMedium:TypeAtrium             -8.667e+00  1.415e+01  -0.612  0.54328   
#> InflHigh:TypeAtrium                1.000e+00  1.415e+01   0.071  0.94397   
#> InflMedium:TypeTerrace            -4.000e+00  1.415e+01  -0.283  0.77871   
#> InflHigh:TypeTerrace               1.667e+00  1.415e+01   0.118  0.90676   
#> InflMedium:ContHigh               -4.000e+00  1.415e+01  -0.283  0.77871   
#> InflHigh:ContHigh                 -8.667e+00  1.415e+01  -0.612  0.54328   
#> TypeApartment:ContHigh             2.200e+01  1.415e+01   1.555  0.12689   
#> TypeAtrium:ContHigh                1.033e+01  1.415e+01   0.730  0.46897   
#> TypeTerrace:ContHigh               2.067e+01  1.415e+01   1.460  0.15098   
#> InflMedium:TypeApartment:ContHigh  2.333e+00  2.001e+01   0.117  0.90769   
#> InflHigh:TypeApartment:ContHigh   -1.200e+01  2.001e+01  -0.600  0.55171   
#> InflMedium:TypeAtrium:ContHigh     3.000e+00  2.001e+01   0.150  0.88150   
#> InflHigh:TypeAtrium:ContHigh       3.667e+00  2.001e+01   0.183  0.85544   
#> InflMedium:TypeTerrace:ContHigh   -8.667e+00  2.001e+01  -0.433  0.66700   
#> InflHigh:TypeTerrace:ContHigh     -1.167e+01  2.001e+01  -0.583  0.56278   
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 12.26 on 46 degrees of freedom
#> Multiple R-squared:  0.6882,	Adjusted R-squared:  0.5187 
#> F-statistic: 4.061 on 25 and 46 DF,  p-value: 1.939e-05
#> 

r.gam <- FitMod(Freq ~ Infl*Type*Cont + Sat, data=housing, fitfn="gamma")
summary(r.gam)
#> 
#> Call:
#> glm(formula = Freq ~ Infl * Type * Cont + Sat, family = "Gamma", 
#>     data = housing)
#> 
#> Coefficients:
#>                                     Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                        4.355e-02  1.143e-02   3.809 0.000412 ***
#> InflMedium                        -1.005e-02  1.433e-02  -0.702 0.486503    
#> InflHigh                           9.653e-03  1.811e-02   0.533 0.596679    
#> TypeApartment                     -1.288e-02  1.386e-02  -0.929 0.357851    
#> TypeAtrium                         5.053e-02  2.761e-02   1.830 0.073681 .  
#> TypeTerrace                        5.354e-02  2.835e-02   1.889 0.065241 .  
#> ContHigh                          -1.378e-17  1.615e-02   0.000 1.000000    
#> Sat.L                             -3.211e-03  3.155e-03  -1.018 0.314213    
#> Sat.Q                             -9.144e-03  4.197e-03  -2.179 0.034492 *  
#> InflMedium:TypeApartment           5.915e-03  1.766e-02   0.335 0.739147    
#> InflHigh:TypeApartment            -8.770e-03  2.134e-02  -0.411 0.683073    
#> InflMedium:TypeAtrium              2.340e-02  4.077e-02   0.574 0.568768    
#> InflHigh:TypeAtrium                3.286e-02  4.794e-02   0.685 0.496500    
#> InflMedium:TypeTerrace            -1.346e-02  3.553e-02  -0.379 0.706642    
#> InflHigh:TypeTerrace               2.393e-02  4.718e-02   0.507 0.614453    
#> InflMedium:ContHigh                4.784e-03  2.086e-02   0.229 0.819606    
#> InflHigh:ContHigh                  4.389e-02  3.364e-02   1.305 0.198497    
#> TypeApartment:ContHigh            -1.121e-02  1.856e-02  -0.604 0.548836    
#> TypeAtrium:ContHigh               -4.583e-02  3.247e-02  -1.412 0.164758    
#> TypeTerrace:ContHigh              -6.393e-02  3.174e-02  -2.014 0.049831 *  
#> InflMedium:TypeApartment:ContHigh -1.765e-03  2.412e-02  -0.073 0.941979    
#> InflHigh:TypeApartment:ContHigh   -3.385e-02  3.662e-02  -0.924 0.360209    
#> InflMedium:TypeAtrium:ContHigh    -1.225e-02  4.752e-02  -0.258 0.797727    
#> InflHigh:TypeAtrium:ContHigh      -5.532e-02  6.091e-02  -0.908 0.368570    
#> InflMedium:TypeTerrace:ContHigh    3.236e-02  4.143e-02   0.781 0.438774    
#> InflHigh:TypeTerrace:ContHigh      1.462e-02  6.501e-02   0.225 0.823077    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for Gamma family taken to be 0.2163103)
#> 
#>     Null deviance: 38.140  on 71  degrees of freedom
#> Residual deviance: 10.671  on 46  degrees of freedom
#> AIC: 538.18
#> 
#> Number of Fisher Scoring iterations: 6
#> 

r.gami <- FitMod(Freq ~ Infl*Type*Cont + Sat, data=housing, fitfn="gamma", link="identity")
summary(r.gami)
#> 
#> Call:
#> glm(formula = Freq ~ Infl * Type * Cont + Sat, family = function () 
#> Gamma(link = identity), data = housing)
#> 
#> Coefficients:
#>                                   Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                        23.3635     5.7553   4.060 0.000189 ***
#> InflMedium                          7.1563     9.5281   0.751 0.456438    
#> InflHigh                           -7.2034     6.9307  -1.039 0.304079    
#> TypeApartment                      12.6892    10.6939   1.187 0.241486    
#> TypeAtrium                        -11.1942     6.3987  -1.749 0.086883 .  
#> TypeTerrace                       -10.9532     6.4273  -1.704 0.095097 .  
#> ContHigh                           -1.6216     7.8406  -0.207 0.837058    
#> Sat.L                               4.9368     1.4858   3.323 0.001755 ** 
#> Sat.Q                               3.2415     1.1421   2.838 0.006731 ** 
#> InflMedium:TypeApartment           -3.4540    16.4721  -0.210 0.834835    
#> InflHigh:TypeApartment              1.4504    13.6446   0.106 0.915810    
#> InflMedium:TypeAtrium              -9.2801    10.1778  -0.912 0.366629    
#> InflHigh:TypeAtrium                 3.5471     7.6884   0.461 0.646714    
#> InflMedium:TypeTerrace             -4.5638    10.5680  -0.432 0.667870    
#> InflHigh:TypeTerrace                2.9450     7.6893   0.383 0.703488    
#> InflMedium:ContHigh                -3.5733    12.5840  -0.284 0.777715    
#> InflHigh:ContHigh                  -7.1152     8.8705  -0.802 0.426610    
#> TypeApartment:ContHigh             22.5209    18.6701   1.206 0.233888    
#> TypeAtrium:ContHigh                11.5520     9.9375   1.162 0.251043    
#> TypeTerrace:ContHigh               22.9328    11.8581   1.934 0.059288 .  
#> InflMedium:TypeApartment:ContHigh   0.9018    27.5339   0.033 0.974015    
#> InflHigh:TypeApartment:ContHigh   -13.2180    21.9476  -0.602 0.549965    
#> InflMedium:TypeAtrium:ContHigh      2.3032    14.8708   0.155 0.877591    
#> InflHigh:TypeAtrium:ContHigh        0.3769    11.2392   0.034 0.973392    
#> InflMedium:TypeTerrace:ContHigh    -8.8401    16.8767  -0.524 0.602931    
#> InflHigh:TypeTerrace:ContHigh     -14.3587    12.7806  -1.123 0.267063    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for Gamma family taken to be 0.1923107)
#> 
#>     Null deviance: 38.1400  on 71  degrees of freedom
#> Residual deviance:  9.0822  on 46  degrees of freedom
#> AIC: 526.3
#> 
#> Number of Fisher Scoring iterations: 14
#> 

old <-options(digits=3)
TMod(r.pois, r.nb, r.log, r.ols, r.gam, r.gami)
#> Waiting for profiling to be done...
#> Waiting for profiling to be done...
#> Waiting for profiling to be done...
#> Waiting for profiling to be done...
#>                                 coef   r.pois         r.nb       r.log    
#> 1                        (Intercept)    3.136 ***    3.140 ***   3.140 ***
#> 2                         InflMedium    0.273 .      0.266       0.260    
#> 3                           InflHigh   -0.205       -0.253      -0.379    
#> 4                      TypeApartment    0.367 *      0.395       0.219    
#> 5                         TypeAtrium   -0.783 ***   -0.776 *    -0.785 *  
#> 6                        TypeTerrace   -0.815 ***   -0.807 *    -0.931 *  
#> 7                           ContHigh    1.409e-15   -0.024      -0.076    
#> 8                              Sat.L    0.116 **     0.169 *     0.200 *  
#> 9                              Sat.Q    0.263 ***    0.257 ***   0.199 *  
#> 10          InflMedium:TypeApartment   -0.118       -0.128       0.049    
#> 11            InflHigh:TypeApartment    0.175        0.151       0.398    
#> 12             InflMedium:TypeAtrium   -0.407       -0.413      -0.400    
#> 13               InflHigh:TypeAtrium   -0.169       -0.130       0.002    
#> 14            InflMedium:TypeTerrace    0.006        0.017       0.143    
#> 15              InflHigh:TypeTerrace   -0.093       -0.070       0.154    
#> 16               InflMedium:ContHigh   -0.140       -0.126      -0.105    
#> 17                 InflHigh:ContHigh   -0.609 *     -0.591      -0.738    
#> 18            TypeApartment:ContHigh    0.503 *      0.522       0.698    
#> 19               TypeAtrium:ContHigh    0.677 *      0.716       0.763    
#> 20              TypeTerrace:ContHigh    1.099 ***    1.150 **    1.114 *  
#> 21 InflMedium:TypeApartment:ContHigh    0.054        0.012      -0.141    
#> 22   InflHigh:TypeApartment:ContHigh    0.146        0.132       0.088    
#> 23    InflMedium:TypeAtrium:ContHigh    0.156        0.141       0.061    
#> 24      InflHigh:TypeAtrium:ContHigh    0.478        0.424       0.503    
#> 25   InflMedium:TypeTerrace:ContHigh   -0.498       -0.517      -0.531    
#> 26     InflHigh:TypeTerrace:ContHigh   -0.447       -0.497      -0.296    
#> 27                               ---                                      
#> 28                         r.squared    -            -           0.751    
#> 29                     adj.r.squared    -            -           0.616    
#> 30                             sigma    -            -           0.463    
#> 31                            logLik -279.213     -240.645     -30.648    
#> 32                           logLik0 -587.313     -290.899       -        
#> 33                                G2  616.201      100.508       -        
#> 34                          deviance    -            -           9.876    
#> 35                               AIC  610.426      535.289     115.296    
#> 36                               BIC  669.619      596.759     176.766    
#> 37                             numdf   26           27          25        
#> 38                             dendf    -            -          46        
#> 39                                 N   72           72          72        
#> 40                            n vars    8            8           8        
#> 41                            n coef   26           26          26        
#> 42                                 F    -            -           5.551    
#> 43                                 p    -            -           0.000    
#> 44                               MAE    6.641        6.608       0.276    
#> 45                              MAPE    0.354        0.345       0.104    
#> 46                               MSE   92.754       97.134       0.137    
#> 47                              RMSE    9.631        9.856       0.370    
#> 48                          McFadden    0.525        0.173       -        
#> 49                       McFaddenAdj    0.480        0.083       -        
#> 50                        Nagelkerke    1.000        0.753       -        
#> 51                          CoxSnell    1.000        0.752       -        
#>       r.ols        r.gam       r.gami    
#> 1    23.333 **     0.044 ***   23.363 ***
#> 2     7.333       -0.010        7.156    
#> 3    -4.333        0.010       -7.203    
#> 4    10.333       -0.013       12.689    
#> 5   -12.667        0.051 .    -11.194 .  
#> 6   -13.000        0.054 .    -10.953 .  
#> 7    -1.680e-14   -0.000       -1.622    
#> 8     2.976       -0.003        4.937 ** 
#> 9     5.835 *     -0.009 *      3.241 ** 
#> 10   -1.667        0.006       -3.454    
#> 11    3.333       -0.009        1.450    
#> 12   -8.667        0.023       -9.280    
#> 13    1.000        0.033        3.547    
#> 14   -4.000       -0.013       -4.564    
#> 15    1.667        0.024        2.945    
#> 16   -4.000        0.005       -3.573    
#> 17   -8.667        0.044       -7.115    
#> 18   22.000       -0.011       22.521    
#> 19   10.333       -0.046       11.552    
#> 20   20.667       -0.064 *     22.933 .  
#> 21    2.333       -0.002        0.902    
#> 22  -12.000       -0.034      -13.218    
#> 23    3.000       -0.012        2.303    
#> 24    3.667       -0.055        0.377    
#> 25   -8.667        0.032       -8.840    
#> 26  -11.667        0.015      -14.359    
#> 27                                       
#> 28    0.688        -            -        
#> 29    0.519        -            -        
#> 30   12.256        -            -        
#> 31 -266.465     -242.088     -236.152    
#> 32    -         -290.204     -290.204    
#> 33    -           96.232      108.104    
#> 34 6909.139        -            -        
#> 35  586.930      538.176      526.304    
#> 36  648.400      599.646      587.774    
#> 37   25           27           27        
#> 38   46            -            -        
#> 39   72           72           72        
#> 40    8            8            8        
#> 41   26           26           26        
#> 42    4.061        -            -        
#> 43    0.000        -            -        
#> 44    7.050        6.883        6.823    
#> 45    0.385        0.371        0.321    
#> 46   95.960       98.993      102.471    
#> 47    9.796        9.950       10.123    
#> 48    -            0.166        0.186    
#> 49    -            0.076        0.097    
#> 50    -            0.737        0.777    
#> 51    -            0.737        0.777    
options(old)


## Ordered Regression

r.polr <- FitMod(Sat ~ Infl + Type + Cont, data=housing, fitfn="polr", weights = Freq)

# multinomial Regression
# r.mult <- FitMod(factor(Sat, ordered=FALSE) ~ Infl + Type + Cont, data=housing,
#                  weights = housing$Freq, fitfn="multinom")


# Regression tree
r.rp <- FitMod(factor(Sat, ordered=FALSE) ~ Infl + Type + Cont, data=housing,
                 weights = housing$Freq, fitfn="rpart")

# compare predictions
d.p <- expand.grid(Infl=levels(housing$Infl), Type=levels(housing$Type), Cont=levels(housing$Cont))
d.p$polr <- predict(r.polr, newdata=d.p)
# ??
# d.p$ols <- factor(round(predict(r.ols, newdata=d.p)^2), labels=levels(housing$Sat))
# d.p$mult <- predict(r.mult, newdata=d.p)
d.p$rp <- predict(r.rp, newdata=d.p, type="class")

d.p
#>      Infl      Type Cont polr   rp
#> 1     Low     Tower  Low  Low High
#> 2  Medium     Tower  Low High High
#> 3    High     Tower  Low High High
#> 4     Low Apartment  Low  Low  Low
#> 5  Medium Apartment  Low  Low High
#> 6    High Apartment  Low High High
#> 7     Low    Atrium  Low  Low High
#> 8  Medium    Atrium  Low High High
#> 9    High    Atrium  Low High High
#> 10    Low   Terrace  Low  Low  Low
#> 11 Medium   Terrace  Low  Low High
#> 12   High   Terrace  Low High High
#> 13    Low     Tower High High High
#> 14 Medium     Tower High High High
#> 15   High     Tower High High High
#> 16    Low Apartment High  Low  Low
#> 17 Medium Apartment High High High
#> 18   High Apartment High High High
#> 19    Low    Atrium High  Low High
#> 20 Medium    Atrium High High High
#> 21   High    Atrium High High High
#> 22    Low   Terrace High  Low  Low
#> 23 Medium   Terrace High  Low High
#> 24   High   Terrace High High High


# Classification with 2 classes  ***************

r.pima <- FitMod(diabetes ~ ., d.pima, fitfn="logit")
r.pima
#> 
#> Call:  glm(formula = diabetes ~ ., family = "binomial", data = d.pima)
#> 
#> Coefficients:
#> (Intercept)     pregnant      glucose     pressure      triceps      insulin  
#>   -8.404696     0.123182     0.035164    -0.013296     0.000619    -0.001192  
#>        mass     pedigree          age  
#>    0.089701     0.945180     0.014869  
#> 
#> Degrees of Freedom: 767 Total (i.e. Null);  759 Residual
#> Null Deviance:	    993.5 
#> Residual Deviance: 723.4 	AIC: 741.4
Conf(r.pima)
#> 
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction pos neg
#>        pos 156  55
#>        neg 112 445
#> 
#>                 Total n : 768
#>                Accuracy : 0.7826
#>                  95% CI : (0.7520, 0.8103)
#>     No Information Rate : 0.6510
#>     P-Value [Acc > NIR] : 1.37e-15
#> 
#>                   Kappa : 0.4966
#>  Mcnemar's Test P-Value : 1.47e-05
#> 
#>             Sensitivity : 0.5821
#>             Specificity : 0.8900
#>          Pos Pred Value : 0.7393
#>          Neg Pred Value : 0.7989
#>              Prevalence : 0.3490
#>          Detection Rate : 0.2747
#>    Detection Prevalence : 0.2031
#>       Balanced Accuracy : 0.7360
#>          F-val Accuracy : 0.6514
#>      Matthews Cor.-Coef : 0.5041
#> 
#>        'Positive' Class : pos
#> 
plot(ROC(r.pima))
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
OddsRatio(r.pima)
#> 
#> Call:
#> glm(formula = diabetes ~ ., family = "binomial", data = d.pima)
#> 
#> Odds Ratios:
#>                or or.lci or.uci  Pr(>|z|)    
#> (Intercept) 0.000  0.000  0.001 < 2.2e-16 ***
#> pregnant    1.131  1.062  1.204  1.23e-04 ***
#> glucose     1.036  1.028  1.043 < 2.2e-16 ***
#> pressure    0.987  0.977  0.997    0.0111 *  
#> triceps     1.001  0.987  1.014    0.9285    
#> insulin     0.999  0.997  1.001    0.1861    
#> mass        1.094  1.062  1.127  2.76e-09 ***
#> pedigree    2.573  1.432  4.625    0.0016 ** 
#> age         1.015  0.997  1.034    0.1112    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Brier Score: 0.153     Nagelkerke R2: 0.408
#> 


# rpart tree
rp.pima <- FitMod(diabetes ~ ., d.pima, fitfn="rpart")
rp.pima
#> n= 768 
#> 
#> node), split, n, loss, yval, (yprob)
#>       * denotes terminal node
#> 
#>   1) root 768 268 neg (0.65104167 0.34895833)  
#>     2) glucose< 127.5 485  94 neg (0.80618557 0.19381443)  
#>       4) age< 28.5 271  23 neg (0.91512915 0.08487085) *
#>       5) age>=28.5 214  71 neg (0.66822430 0.33177570)  
#>        10) mass< 26.35 41   2 neg (0.95121951 0.04878049) *
#>        11) mass>=26.35 173  69 neg (0.60115607 0.39884393)  
#>          22) glucose< 99.5 55  10 neg (0.81818182 0.18181818) *
#>          23) glucose>=99.5 118  59 neg (0.50000000 0.50000000)  
#>            46) pedigree< 0.561 84  34 neg (0.59523810 0.40476190)  
#>              92) pedigree< 0.2 21   4 neg (0.80952381 0.19047619) *
#>              93) pedigree>=0.2 63  30 neg (0.52380952 0.47619048)  
#>               186) pregnant>=1.5 52  21 neg (0.59615385 0.40384615)  
#>                 372) pressure>=67 40  12 neg (0.70000000 0.30000000) *
#>                 373) pressure< 67 12   3 pos (0.25000000 0.75000000) *
#>               187) pregnant< 1.5 11   2 pos (0.18181818 0.81818182) *
#>            47) pedigree>=0.561 34   9 pos (0.26470588 0.73529412) *
#>     3) glucose>=127.5 283 109 pos (0.38515901 0.61484099)  
#>       6) mass< 29.95 76  24 neg (0.68421053 0.31578947)  
#>        12) glucose< 145.5 41   6 neg (0.85365854 0.14634146) *
#>        13) glucose>=145.5 35  17 pos (0.48571429 0.51428571)  
#>          26) insulin< 14.5 21   8 neg (0.61904762 0.38095238) *
#>          27) insulin>=14.5 14   4 pos (0.28571429 0.71428571) *
#>       7) mass>=29.95 207  57 pos (0.27536232 0.72463768)  
#>        14) glucose< 157.5 115  45 pos (0.39130435 0.60869565)  
#>          28) age< 30.5 50  23 neg (0.54000000 0.46000000)  
#>            56) pressure>=61 40  13 neg (0.67500000 0.32500000)  
#>             112) mass< 41.8 31   7 neg (0.77419355 0.22580645) *
#>             113) mass>=41.8 9   3 pos (0.33333333 0.66666667) *
#>            57) pressure< 61 10   0 pos (0.00000000 1.00000000) *
#>          29) age>=30.5 65  18 pos (0.27692308 0.72307692) *
#>        15) glucose>=157.5 92  12 pos (0.13043478 0.86956522) *
Conf(rp.pima)
#> 
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction neg pos
#>        neg 449  72
#>        pos  51 196
#> 
#>                 Total n : 768
#>                Accuracy : 0.8398
#>                  95% CI : (0.8122, 0.8641)
#>     No Information Rate : 0.6510
#>     P-Value [Acc > NIR] : < 2.2e-16
#> 
#>                   Kappa : 0.6410
#>  Mcnemar's Test P-Value : 0.0713
#> 
#>             Sensitivity : 0.8980
#>             Specificity : 0.7313
#>          Pos Pred Value : 0.8618
#>          Neg Pred Value : 0.7935
#>              Prevalence : 0.6510
#>          Detection Rate : 0.6784
#>    Detection Prevalence : 0.5846
#>       Balanced Accuracy : 0.8147
#>          F-val Accuracy : 0.8795
#>      Matthews Cor.-Coef : 0.6422
#> 
#>        'Positive' Class : neg
#> 
lines(ROC(rp.pima), col=hblue)
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases

# to be improved
plot(rp.pima, col=SetAlpha(c("blue","red"), 0.4), cex=0.7)


# Random Forest
rf.pima <- FitMod(diabetes ~ ., d.pima, method="class", fitfn="randomForest")
rf.pima
#> 
#> Call:
#>  randomForest(formula = diabetes ~ ., data = d.pima, method = "class",      na.action = function (object, ...)      UseMethod("na.omit")) 
#>                Type of random forest: classification
#>                      Number of trees: 500
#> No. of variables tried at each split: 2
#> 
#>         OOB estimate of  error rate: 23.05%
#> Confusion matrix:
#>     neg pos class.error
#> neg 430  70   0.1400000
#> pos 107 161   0.3992537
Conf(rf.pima)
#> 
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction neg pos
#>        neg 430 107
#>        pos  70 161
#> 
#>                 Total n : 768
#>                Accuracy : 0.7695
#>                  95% CI : (0.7384, 0.7979)
#>     No Information Rate : 0.6510
#>     P-Value [Acc > NIR] : 7.16e-13
#> 
#>                   Kappa : 0.4760
#>  Mcnemar's Test P-Value : 0.0068
#> 
#>             Sensitivity : 0.8600
#>             Specificity : 0.6007
#>          Pos Pred Value : 0.8007
#>          Neg Pred Value : 0.6970
#>              Prevalence : 0.6510
#>          Detection Rate : 0.6992
#>    Detection Prevalence : 0.5599
#>       Balanced Accuracy : 0.7304
#>          F-val Accuracy : 0.8293
#>      Matthews Cor.-Coef : 0.4789
#> 
#>        'Positive' Class : neg
#> 
lines(ROC(r.pima), col=hred)
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases




# more models to compare

d.pim <- SplitTrainTest(d.pima, p = 0.2)
mdiab <- formula(diabetes ~ pregnant + glucose + pressure + triceps
                            + insulin + mass + pedigree + age)

r.glm <- FitMod(mdiab, data=d.pim$train, fitfn="logit")
r.rp <- FitMod(mdiab, data=d.pim$train, fitfn="rpart")
r.rf <- FitMod(mdiab, data=d.pim$train, fitfn="randomForest")
r.svm <- FitMod(mdiab, data=d.pim$train, fitfn="svm")
r.c5 <- FitMod(mdiab, data=d.pim$train, fitfn="C5.0")
r.nn <- FitMod(mdiab, data=d.pim$train, fitfn="nnet")
r.nb <- FitMod(mdiab, data=d.pim$train, fitfn="naive_bayes")
r.lda <- FitMod(mdiab, data=d.pim$train, fitfn="lda")
r.qda <- FitMod(mdiab, data=d.pim$train, fitfn="qda")
r.lb <- FitMod(mdiab, data=d.pim$train, fitfn="lb")

mods <- list(glm=r.glm, rp=r.rp, rf=r.rf, svm=r.svm, c5=r.c5
             , nn=r.nn, nb=r.nb, lda=r.lda, qda=r.qda, lb=r.lb)

# insight in the Regression tree
plot(r.rp, box.palette = as.list(Pal("Helsana", alpha = 0.5)))


# Insample accuracy ...
TModC(mods, ord="auc")
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Error in eval(object$call$data): object 'd.pim' not found
# ... is substantially different from the out-of-bag:
TModC(mods, newdata=d.pim$test, reference=d.pim$test$diabetes, ord="bs")
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Warning: predict.naive_bayes(): more features in the newdata are provided as there are probability tables in the object. Calculation is performed based on features to be found in the tables.
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Setting levels: control = neg, case = pos
#> Setting direction: controls < cases
#> Error in eval(x$call$formula): object 'mdiab' not found
# C5 and SVM turn out to be show-offs! They overfit quite ordinary
# whereas randomforest and logit keep their promises. ...

sapply(mods, function(z) VarImp(z))
#> Warning: no VarImp definition found for x (class FitMod, naive_bayes)
#> Warning: no VarImp definition found for x (class FitMod, lda)
#> Warning: no VarImp definition found for x (class FitMod, qda)
#> Warning: no VarImp definition found for x (class FitMod, LogitBoost)
#> $glm
#>          Overall
#> pregnant   3.210
#> glucose    8.176
#> pressure   2.437
#> triceps    0.056
#> insulin    1.642
#> mass       5.010
#> pedigree   2.970
#> age        1.687
#> 
#> $rp
#>          varimp
#> glucose  57.127
#> mass     17.811
#> pressure 13.715
#> age      13.531
#> insulin  12.599
#> pregnant  2.970
#> triceps   2.931
#> pedigree  0.866
#> 
#> $rf
#>          varimp
#> glucose  68.585
#> mass     44.600
#> age      38.308
#> pedigree 35.697
#> pressure 25.106
#> pregnant 24.142
#> triceps  20.404
#> insulin  18.435
#> 
#> $svm
#>      x
#> 1 <NA>
#> 
#> $c5
#>          Overall
#> glucose  100.000
#> mass     100.000
#> age       66.180
#> triceps   20.000
#> pregnant  18.860
#> pedigree  15.610
#> pressure   7.320
#> insulin    0.000
#> 
#> $nn
#>          Overall
#> triceps   17.802
#> pedigree  16.997
#> glucose   15.842
#> pressure  14.737
#> pregnant  12.023
#> age        8.704
#> insulin    7.800
#> mass       6.095
#> 
#> $nb
#>      x
#> 1 <NA>
#> 
#> $lda
#>      x
#> 1 <NA>
#> 
#> $qda
#>      x
#> 1 <NA>
#> 
#> $lb
#>      x
#> 1 <NA>
#> 


# Multinomial classification problem with n classes  ***************

d.gl <- SplitTrainTest(d.glass, p = 0.2)
mglass <- formula(Type ~ RI + Na + Mg + Al + Si + K + Ca + Ba + Fe)

# *** raises an unclear error in CRAN-Debian tests *** ??
# r.mult <- FitMod(mglass, data=d.gl$train, maxit=600, fitfn="multinom")
r.rp <- FitMod(mglass, data=d.gl$train, fitfn="rpart")
r.rf <- FitMod(mglass, data=d.gl$train, fitfn="randomForest")
r.svm <- FitMod(mglass, data=d.gl$train, fitfn="svm")
r.c5 <- FitMod(mglass, data=d.gl$train, fitfn="C5.0")
r.nn <- FitMod(mglass, data=d.gl$train, fitfn="nnet")
r.nbay <- FitMod(mglass, data=d.gl$train, fitfn="naive_bayes")
r.lda <- FitMod(mglass, data=d.gl$train, fitfn="lda")
# r.qda <- FitMod(mglass, data=d.glass, fitfn="qda")
r.lb <- FitMod(mglass, data=d.gl$train, fitfn="lb")

mods <- list(rp=r.rp, rf=r.rf, svm=r.svm, c5=r.c5,
             nn=r.nn, nbay=r.nbay, lda=r.lda, lb=r.lb)

# confusion matrix and other quality measures can be calculated with Conf()
Conf(r.rf)
#> 
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction  1  2  3  5  6  7
#>          1 51 10  6  0  0  1
#>          2  6 52  2  3  1  3
#>          3  1  0  5  0  0  0
#>          5  0  1  0  6  0  0
#>          6  0  0  0  0  6  0
#>          7  0  1  0  1  0 16
#> 
#> Overall Statistics
#> 
#>                 Total n : 172
#>                Accuracy : 0.7907
#>                  95% CI : (0.7239, 0.8448)
#>     No Information Rate : 0.3721
#>     P-Value [Acc > NIR] : < 2.2e-16
#> 
#>                   Kappa : 0.7023
#>  Mcnemar's Test P-Value : NA
#> 
#> 
#> Statistics by Class:
#> 
#>                           1      2      3      5      6      7
#> Sensitivity          0.8793 0.8125 0.3846 0.6000 0.8571 0.8000
#> Specificity          0.8509 0.8611 0.9937 0.9938 1.0000 0.9868
#> Pos Pred Value       0.7500 0.7761 0.8333 0.8571 1.0000 0.8889
#> Neg Pred Value       0.9327 0.8857 0.9518 0.9758 0.9940 0.9740
#> Prevalence           0.3372 0.3721 0.0756 0.0581 0.0407 0.1163
#> Detection Rate       0.3953 0.3895 0.0349 0.0407 0.0349 0.1047
#> Detection Prevalence 0.2965 0.3023 0.0291 0.0349 0.0349 0.0930
#> Balanced Accuracy    0.8651 0.8368 0.6892 0.7969 0.9286 0.8934
#> F-val Accuracy       0.8095 0.7939 0.5263 0.7059 0.9231 0.8421
#> Matthews Cor.-Coef   0.7060 0.6677 0.5450 0.7033 0.9230 0.8240
#> 

# we only extract the general accuracy
sapply(lapply(mods, function(z) Conf(z)), "[[", "acc")
#> Error in eval(object$call$data): object 'd.gl' not found

# let's compare r.mult with a model without RI as predictor
# Conf(r.mult)
# Conf(update(r.mult, . ~ . -RI))