Construire un modèle parsnip

Introduction

Notre objectif ici est de construire un nouveau modèle {parsnip} à partir d’une implémentation existante dans un package R.

Intérêts : éliminer des dépendances, du code dupliqué, pouvoir profiter du cadre {tidymodels}.

Ressources

• https://www.tidymodels.org/learn/develop/models/

• Tuto tidymodels - Rencontres R 2023

Un modèle {parsnip}

Il se définit par

• un mode (régression linéaire, logistique, …),

• un type (régression, classification),

• un moteur de calcul engine (lm, stan, …).

Qunad on ajoute un modèle, on doit donc spécifier son mode et son engine. On peut aussi ajouter un nouveau modèle différent d’un pré-existant seulement par son type (ie même combinaison engine et mode).

Il est possible de voir l’ensemble des modèles déjà disponibles ici.

Intégration de la régression PLN

library(PLNmodels)
library(tidymodels)
library(poissonreg)

Données

Nous utiliserons le jeu de données de la vignette de PLN.

data(trichoptera)

Remarque

Il existe une fonction prepare_data dans {parsnip} également.

Différentes étapes pour intégrer PLNmodels::PLN dans parsnip

Etape 1 : Spécification du modèle et de ses arguments

On va ajouter la fonction PLNmodels::PLN à la liste des fonctions utilisables pour faire de la régression de Poisson.

show_model_info("poisson_reg")

Information for `poisson_reg`
 modes: unknown, regression 

 engines: 
   regression: glm¹, glmnet¹, hurdle¹, stan¹, zeroinfl¹

¹The model can use case weights.

 arguments: 
   glmnet: 
      penalty --> lambda
      mixture --> alpha

 fit modules:
     engine       mode
        glm regression
     hurdle regression
   zeroinfl regression
     glmnet regression
       stan regression

 prediction modules:
         mode   engine                          methods
   regression      glm                     numeric, raw
   regression   glmnet                     numeric, raw
   regression   hurdle                     numeric, raw
   regression     stan conf_int, numeric, pred_int, raw
   regression zeroinfl                     numeric, raw

set_model_mode(model = "poisson_reg", mode = "regression")
set_model_engine(
  "poisson_reg", 
  mode = "regression", 
  eng = "PLN"
)
set_dependency("poisson_reg", eng = "PLN", pkg = "PLNmodels")

On peut vérifier ce qu’on a ajouté.

show_model_info("poisson_reg")

Information for `poisson_reg`
 modes: unknown, regression 

 engines: 
   regression: glm¹, glmnet¹, hurdle¹, PLNNA, stan¹, zeroinfl¹

¹The model can use case weights.

 arguments: 
   glmnet: 
      penalty --> lambda
      mixture --> alpha

 fit modules:
     engine       mode
        glm regression
     hurdle regression
   zeroinfl regression
     glmnet regression
       stan regression

 prediction modules:
         mode   engine                          methods
   regression      glm                     numeric, raw
   regression   glmnet                     numeric, raw
   regression   hurdle                     numeric, raw
   regression     stan conf_int, numeric, pred_int, raw
   regression zeroinfl                     numeric, raw

Si notre modèle a des arguments supplémentaires que ceux existants, on peut les ajouter.

Etape 2 : Créer la fonction principale associée au modèle le cas échéant

Ici comme nous avons seulement ajouté un moteur à un modèle préexistant ce n’est pas nécessaire.

todo: regarder comment ajouter l’information possible d’utilisation de poids

Etape 3 : Préciser le module d’ajustement (fit)

• L’argument interface, peut prendre les valeurs formula, data ou matrices.

• protect est une liste optionnelle d’arguments qui ne devraient pas être modifié par l’utilisateur.

• func est un vecteur précisant le package et la fonction qui sera appelé pour l’ajustement.

• defaults est une liste optionnelle d’arguments que l’utilisateur peut modifier mais pour laquelle on peut spécifier des valeurs par défaut.

set_fit(
  model = "poisson_reg",
  eng = "PLN",
  mode = "regression",
  value = list(
    interface = "formula",
    protect = c("formula", "data"),
    func = c(pkg = "PLNmodels", fun = "PLN"),
    defaults = list(control = PLN_param(), weights = NULL, subset = NULL)
  )
)

On peut ajouter des traitements par défaut pour les variables explicatives, tels que calculs de variables indicatrices, calcul d’une constante etc…

set_encoding(
  model = "poisson_reg",
  eng = "PLN",
  mode = "regression",
  options = list(
    predictor_indicators = "traditional",
    compute_intercept = TRUE,
    remove_intercept = TRUE,
    allow_sparse_x = FALSE
  )
)

Etape 4 : Ajouter un module pour la prédiction (predict)

set_pred(
  model = "poisson_reg",
  eng = "PLN",
  mode = "regression",
  type = "numeric",
  value = list(
    pre = NULL,
    post = NULL,
    func = c(fun = "predict"),
    args =
      list(
        object = expr(object$fit),
        newdata = expr(new_data),
        type = "response"
      )
  )
)

Application sur les données trichoptera

PLN dans PLNmodels

prep_trichoptera <- PLNmodels::prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1 , data = prep_trichoptera)

 Initialization...
 Adjusting a full covariance PLN model with nlopt optimizer
 Post-treatments...
 DONE!

myPLN

A multivariate Poisson Lognormal fit with full covariance model.
==================================================================
 nb_param    loglik       BIC       ICL
      170 -1129.909 -1460.714 -2230.658
==================================================================
* Useful fields
    $model_par, $latent, $latent_pos, $var_par, $optim_par
    $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), vcov(), fitted()
    predict(), predict_cond(), standard_error()

PLN dans parsnip

resPLN <- poisson_reg() %>% 
  set_engine("PLN") %>% 
  fit(Abundance ~ 1 , data = prep_trichoptera)

 Initialization...
 Adjusting a full covariance PLN model with nlopt optimizer
 Post-treatments...
 DONE!

resPLN

parsnip model object

A multivariate Poisson Lognormal fit with full covariance model.
==================================================================
 nb_param    loglik       BIC       ICL
      170 -1129.909 -1460.714 -2230.658
==================================================================
* Useful fields
    $model_par, $latent, $latent_pos, $var_par, $optim_par
    $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), vcov(), fitted()
    predict(), predict_cond(), standard_error()

summary(resPLN)

             Length Class       Mode       
lvl           0     -none-      NULL       
spec          8     poisson_reg list       
fit          32     PLNfit      environment
preproc       1     -none-      list       
elapsed       1     -none-      list       
censor_probs  0     -none-      list       

resPLN$spec

Poisson Regression Model Specification (regression)

Computational engine: PLN 

Model fit template:
PLNmodels::PLN(formula = missing_arg(), data = missing_arg(), 
    control = list(backend = "nlopt", trace = 1, covariance = "full", 
        Omega = NULL, config_post = list(jackknife = FALSE, bootstrap = 0L, 
            rsquared = TRUE, variational_var = FALSE, sandwich_var = FALSE, 
            trace = 1), config_optim = list(algorithm = "CCSAQ", 
            maxeval = 10000, ftol_rel = 1e-08, xtol_rel = 1e-06, 
            ftol_abs = 0, xtol_abs = 0, maxtime = -1, trace = 1), 
        inception = NULL), weights = NULL, subset = NULL)

# Pour recuperer les coefficients
coef(resPLN$fit)

                  Che       Hyc      Hym      Hys      Psy        Aga       Glo
(Intercept) -2.833822 -4.049306 1.051399 -2.17906 3.291263 -0.3560404 -2.217591
                  Ath      Cea       Ced        Set        All       Han
(Intercept) -1.661659 -3.36122 0.3593471 -0.6285713 -0.8500234 -1.351443
                  Hfo        Hsp       Hve       Sta
(Intercept) -2.190792 -0.3070371 -2.758038 0.9621415

all.equal(coef(myPLN), coef(resPLN$fit))

[1] TRUE

# Pour faire de la prediction
pred_pln <- predict(myPLN, newdata  = prep_trichoptera, type = "response")
pred_pln %>% head(n=3)

         Che       Hyc      Hym       Hys      Psy      Aga      Glo       Ath
1 0.06157233 0.0588531 3.813198 0.1225141 89.53902 2.888241 0.251068 0.2913597
2 0.06157233 0.0588531 3.813198 0.1225141 89.53902 2.888241 0.251068 0.2913597
3 0.06157233 0.0588531 3.813198 0.1225141 89.53902 2.888241 0.251068 0.2913597
        Cea     Ced      Set      All      Han      Hfo      Hsp       Hve
1 0.1053727 2.54603 3.594756 1.074918 11.25454 1.024682 5.592037 0.1974294
2 0.1053727 2.54603 3.594756 1.074918 11.25454 1.024682 5.592037 0.1974294
3 0.1053727 2.54603 3.594756 1.074918 11.25454 1.024682 5.592037 0.1974294
       Sta
1 7.385873
2 7.385873
3 7.385873

# On peut appliquer les methodes 
#pred_parsnipfit <- predict(resPLN$fit, newdata  = prep_trichoptera)
pred_parsnip <- predict(resPLN, new_data  = trichoptera$Abundance)
pred_parsnip %>% head(n=3)

# A tibble: 3 × 17
  .pred_Che .pred_Hyc .pred_Hym .pred_Hys .pred_Psy .pred_Aga .pred_Glo
      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
1    0.0616    0.0589      3.81     0.123      89.5      2.89     0.251
2    0.0616    0.0589      3.81     0.123      89.5      2.89     0.251
3    0.0616    0.0589      3.81     0.123      89.5      2.89     0.251
# ℹ 10 more variables: .pred_Ath <dbl>, .pred_Cea <dbl>, .pred_Ced <dbl>,
#   .pred_Set <dbl>, .pred_All <dbl>, .pred_Han <dbl>, .pred_Hfo <dbl>,
#   .pred_Hsp <dbl>, .pred_Hve <dbl>, .pred_Sta <dbl>

Test d’un workflow

# separation du jeu de donnees en test et apprentissage
set.seed(123)
tri_split <- initial_split(data = prep_trichoptera, prop = 0.9)
#tri_train <- training(tri_split)
#tri_test <- testing(tri_split)
p_recipe <- 
  recipe(Abundance ~ 1 + Temperature, data = training(tri_split))

pln_model <- 
  poisson_reg()%>%
  set_engine("PLN") %>%
  set_mode("regression")

pln_workflow <- workflow() %>%
  add_recipe(p_recipe) %>%
  add_model(pln_model)

fitted_workflow <-  pln_workflow %>%
  fit(training(tri_split)) 

# Predicition sur le jeu d'entrainement
test_pred <- fitted_workflow %>% predict(new_data = testing(tri_split))

Pour la suite du workflow, calcul de performance sur un rééchantillonnage etc., ce serait possible mais nécessite un peu de code. Il faudrait ici que la sortie de prédiction soit en format plus ‘tidy’ ou réécrire le calcul de métrique sur des matrices.

Autre jeu de données, exemple en régression univariée

p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
p <- within(p, {
  prog <- factor(prog, levels=1:3, labels=c("General", "Academic", 
                                            "Vocational"))
  id <- factor(id)
})

On essaie maintenant en utilisant {parsnip}.

glm de base

summary(m1 <- glm(num_awards ~ prog + math, family="poisson", data=p))

Call:
glm(formula = num_awards ~ prog + math, family = "poisson", data = p)

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -5.24712    0.65845  -7.969 1.60e-15 ***
progAcademic    1.08386    0.35825   3.025  0.00248 ** 
progVocational  0.36981    0.44107   0.838  0.40179    
math            0.07015    0.01060   6.619 3.63e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 287.67  on 199  degrees of freedom
Residual deviance: 189.45  on 196  degrees of freedom
AIC: 373.5

Number of Fisher Scoring iterations: 6

glm dans parsnip

poisson_reg() %>% 
     set_engine("glm") %>%
     set_mode("regression") %>%
     fit(num_awards ~ prog + math, data = p)  %>%
     predict(p)

# A tibble: 200 × 1
    .pred
    <dbl>
 1 0.135 
 2 0.0934
 3 0.167 
 4 0.145 
 5 0.126 
 6 0.100 
 7 0.192 
 8 0.126 
 9 0.0771
10 0.192 
# ℹ 190 more rows

pln

mod_pln <- PLN(as.matrix(num_awards) ~ prog + math, data = p)  

 Initialization...
 Adjusting a full covariance PLN model with nlopt optimizer
 Post-treatments...
 DONE!

mod_pln %>%
     predict(p, type = "response") %>%
     head(n = 10)

            Y
1  0.13408541
2  0.09317098
3  0.16585222
4  0.14393335
5  0.12491126
6  0.10001395
7  0.19110900
8  0.12491126
9  0.07605752
10 0.19110900

PLN dans parsnip

poisson_reg() %>% 
     set_engine("PLN") %>%
     set_mode("regression") %>%
     fit(as.matrix(num_awards) ~ prog + math, data = p) %>% 
     predict(p)

 Initialization...
 Adjusting a full covariance PLN model with nlopt optimizer
 Post-treatments...
 DONE!

# A tibble: 200 × 1
   .pred_Y
     <dbl>
 1  0.134 
 2  0.0932
 3  0.166 
 4  0.144 
 5  0.125 
 6  0.100 
 7  0.191 
 8  0.125 
 9  0.0761
10  0.191 
# ℹ 190 more rows

On peut vérifier la commande associée au modèle programmé.

poisson_reg() %>% translate(engine = "PLN")

Poisson Regression Model Specification (regression)

Computational engine: PLN 

Model fit template:
PLNmodels::PLN(formula = missing_arg(), data = missing_arg(), 
    control = list(backend = "nlopt", trace = 1, covariance = "full", 
        Omega = NULL, config_post = list(jackknife = FALSE, bootstrap = 0L, 
            rsquared = TRUE, variational_var = FALSE, sandwich_var = FALSE, 
            trace = 1), config_optim = list(algorithm = "CCSAQ", 
            maxeval = 10000, ftol_rel = 1e-08, xtol_rel = 1e-06, 
            ftol_abs = 0, xtol_abs = 0, maxtime = -1, trace = 1), 
        inception = NULL), weights = NULL, subset = NULL)

Pour aller plus loin

Dans le cadre tidymodels, tout est assez modulaire et personnalisable : on peut ajouter de nouvelles recettes de pré-traitement, personnaliser son rééchantillonnage, le calcul de métrique etc.

Un peu de documentation sur le site de {tidymodels] dans la section Develop custom modeling tools.