Les paramètres suivants doivent être définis dans R ou RStudio afin de reproduire les analyses présentées dans ce document. Les packages ggplot2, ggfortify, hrbrthemes, directlabels, cowplot, texreg, Hmisc, lme4, geepack et reshape2 ne font pas partie de la distribution R et doivent être installés si nécessaire :

install.packages(c("ggplot2", "ggfortify", "hrbrthemes", "directlabels",
                     "cowplot", "texreg", "Hmisc", "lme4", "geepack", "reshape2")

Les dépendances de ces packages seront installées automatiquement. On supposera également que les instructions R sont exécutées dans un répertoire de travail avec les fichiers de données accessibles dans un sous-répertoire data/. Si ce n’est pas le cas, il suffit de créer un sous-répertoire data/ et d’y enregistrer les fichiers de données, ou de redéfinir les chemins d’accès dans les instructions de lecture des fichiers de données ci-après.

library(ggplot2)
library(ggfortify)
library(hrbrthemes)
library(directlabels)
library(cowplot)
library(texreg)
library(Hmisc)
library(survival)
library(lme4)
library(geepack)
library(reshape2)
options(digits = 6, show.signif.stars = FALSE)
theme_set(theme_ipsum(base_size = 11))

Modèles à effets aléatoires

Chargement des données

Les données analysées ci-après proviennent d’un plan d’expérience dans lequel on s’intéresse à l’excès de graisses dans les selles causé par un déficit en enzymes digestives au niveau de l’intestin. Les mêmes sujets ont été soumis à différents traitements (suppléments en enzyme pancréatique). Il y a donc plusieurs observations par sujet. Les données sont disponibles au format RData dans le fichier fat.rda :

load("data/fat.rda")
str(fat)

## 'data.frame':    24 obs. of  3 variables:
##  $ fecfat  : num  44.5 33 19.1 9.4 71.3 51.2 7.3 21 5 4.6 ...
##  $ pilltype: Factor w/ 4 levels "None","Tablet",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ subject : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 6 1 2 3 4 ...

Description des variables

Les moyennes et écarts-type pour chaque groupe de traitement sont reportées ci-dessous à l’aide de Hmisc :

summary(fecfat ~ pilltype, data = fat, fun = smean.sd)

## fecfat     N= 24 
## 
## +--------+-------+--+-------+-------+
## |        |       |N |Mean   |SD     |
## +--------+-------+--+-------+-------+
## |pilltype|None   | 6|38.0833|22.4745|
## |        |Tablet | 6|16.5333|13.3209|
## |        |Capsule| 6|17.4167|12.9375|
## |        |Coated | 6|31.0667|24.2641|
## +--------+-------+--+-------+-------+
## |Overall |       |24|25.7750|20.0021|
## +--------+-------+--+-------+-------+

Une représentation graphique des données individuelles et moyennes est données ci-après.

d <- aggregate(fecfat ~ pilltype, data = fat, mean)
p <- ggplot(data = fat, aes(x = reorder(pilltype, fecfat), y = fecfat)) +
       geom_line(aes(group = subject), color = grey(.3)) +
       geom_line(data = d, aes(x = pilltype, y = fecfat, group = 1),
                 color = "lightcoral", size = 1.2) +
       labs(x = NULL, y = "Fecal Fat")
p

Modélisation

Voici différents modèles permettant de décomposer la variance totale, dont un modèle d’ANOVA à mesures répétées (modèle m3 ci-dessous) supposant la symétrie composée de la matrice de variance-covariance :

m1 <- aov(fecfat ~ pilltype, data = fat)
m2 <- aov(fecfat ~ pilltype + subject, data = fat)
m3 <- aov(fecfat ~ pilltype + Error(subject), data = fat)
summary(m3)

## 
## Error: subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  5   5588    1118               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## pilltype   3   2009     670    6.26 0.0057
## Residuals 15   1605     107

Et voici une approche par modèle de régression à effets aléatoires, ici un intercept aléatoire pour chaque sujet :

m <- lmer(fecfat ~ pilltype + (1 | subject), data = fat)
summary(m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: fecfat ~ pilltype + (1 | subject)
##    Data: fat
## 
## REML criterion at convergence: 169.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2101 -0.6151 -0.0027  0.4571  1.7256 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 253      15.9    
##  Residual             107      10.3    
## Number of obs: 24, groups:  subject, 6
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)        38.08       7.74    4.92
## pilltypeTablet    -21.55       5.97   -3.61
## pilltypeCapsule   -20.67       5.97   -3.46
## pilltypeCoated     -7.02       5.97   -1.17
## 
## Correlation of Fixed Effects:
##             (Intr) plltyT plltypCp
## pilltypTblt -0.386                
## pilltypCpsl -0.386  0.500         
## pilltypeCtd -0.386  0.500  0.500

confint(m)

## Computing profile confidence intervals ...

##                     2.5 %   97.5 %
## .sig01            8.10894  30.2338
## .sigma            7.03296  13.6222
## (Intercept)      22.52784  53.6388
## pilltypeTablet  -32.83157 -10.2684
## pilltypeCapsule -31.94823  -9.3851
## pilltypeCoated  -18.29823   4.2649

Visuellement, le modèle à intercept aléatoire se traduit par des profils individuels en tous points identiques à l’exception du niveau moyen qui varie d’un individu à l’autre :

yhat <- predict(m)
p <- ggplot(data = fat, aes(x = reorder(pilltype, fecfat), y = yhat)) +
       geom_line(aes(group = subject), color = grey(.3)) +
       geom_line(data = d, aes(x = pilltype, y = fecfat, group = 1),
                 color = "lightcoral", size = 1.2) +
       labs(x = NULL, y = "Predicted Fecal Fat")
p

Cas des données longitudinales

Chargement des données

Les données pour cette illustration portent sur un essai clinique randomisé comprenant deux bras de traitement et visant à étudier l’effet de l’administration d’ibuprofène par voie intraveineuse sur la mortalité de patients en état septique sévère. Les données sont disponibles dans le fichier Stata sepsis :

sepsis <- foreign::read.dta("data/sepsis.dta")
str(sepsis)

## 'data.frame':    455 obs. of  23 variables:
##  $ id      : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ treat   : Factor w/ 2 levels "Placebo","Ibuprofen": 1 2 1 2 1 2 2 1 2 1 ...
##  $ race    : Factor w/ 3 levels "White","Black",..: 1 2 2 1 1 1 1 1 1 1 ...
##  $ apache  : num  27 14 33 3 5 13 34 11 25 20 ...
##  $ o2del   : num  539 NA 551 1376 NA ...
##  $ fate    : Factor w/ 2 levels "Alive","Dead": 2 1 2 1 1 1 2 1 1 1 ...
##  $ followup: num  50 720 33 720 720 720 45 720 720 720 ...
##  $ temp0   : num  95.4 101.6 101 101 101.4 ...
##  $ temp1   : num  93.9 99 98.9 100.2 101.4 ...
##  $ temp2   : num  93.6 99.8 97.2 100 101.8 ...
##  $ temp3   : num  93 99.8 NA 99.5 101.4 ...
##  $ temp4   : num  90.3 99.6 NA 99.7 101.4 ...
##  $ temp5   : num  93.2 100.2 94.8 99.1 101 ...
##  $ temp6   : num  95 100.2 95.5 99.3 102.6 ...
##  $ temp7   : num  96.2 99 NA 99.6 100.8 ...
##  $ temp8   : num  98.8 99.9 97.2 98 100.6 ...
##  $ temp9   : num  97.5 99.8 99.2 97.8 99.8 ...
##  $ temp10  : num  97.8 99.6 NA 97.6 99.6 ...
##  $ temp11  : num  99.1 99.4 NA 98.4 99.6 ...
##  $ temp12  : num  95.6 101.2 NA 98 99.8 ...
##  $ temp13  : num  NA 101 NA 98.4 100.8 ...
##  $ temp14  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ temp15  : num  NA 100.6 NA 98.8 100 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr " 2 Sep 2008 10:17"
##  - attr(*, "formats")= chr  "%9.0g" "%9.0g" "%9.0g" "%9.0g" ...
##  - attr(*, "types")= int  254 254 254 254 254 254 254 254 254 254 ...
##  - attr(*, "val.labels")= chr  "" "treatmnt" "race" "" ...
##  - attr(*, "var.labels")= chr  "Patient ID" "Treatment" "Race" "Baseline APACHE Score" ...
##  - attr(*, "version")= int 12
##  - attr(*, "label.table")=List of 3
##   ..$ treatmnt: Named int  0 1
##   .. ..- attr(*, "names")= chr  "Placebo" "Ibuprofen"
##   ..$ fate    : Named int  0 1
##   .. ..- attr(*, "names")= chr  "Alive" "Dead"
##   ..$ race    : Named int  0 1 2
##   .. ..- attr(*, "names")= chr  "White" "Black" "Other"

Voici quelques pré-traitements utiles :

sepsis$treat <- as.factor(sepsis$treat)
sepsis$id <- factor(sepsis$id)
sepsis.long <- melt(sepsis, measure.vars = 8:14, id.vars = 1:2)
celsius2fahr <- function(x) (x-32) / 1.8              ## °F -> °C
sepsis.long$value <- celsius2fahr(sepsis.long$value)

Description des variables

On ne s’intéresse qu’aux variables treat (groupe de traitement), race (ethnicité), fate (statut vital) et les différentes mesures de température (temp*). Il est possible d’aborder ce jeu de données sous un angle d’analyse de survie, mais on va principalement regarder l’évolution de la température après la prise en charge entre les deux groupes de patients. Voici quelques statistiques descriptives obtenues avec le package Hmisc :

s <- summary(treat ~ fate + race + apache + temp0 + temp1, data = sepsis,
             method = "reverse", overall = TRUE)
print(s, digits = 3)

## 
## 
## Descriptive Statistics by treat
## 
## +------------+---+-----------------+-----------------+-----------------+
## |            |N  |Placebo          |Ibuprofen        |Combined         |
## |            |   |(N=231)          |(N=224)          |(N=455)          |
## +------------+---+-----------------+-----------------+-----------------+
## |fate : Dead |455|    40%  ( 92)   |    38%  ( 84)   |    39%  (176)   |
## +------------+---+-----------------+-----------------+-----------------+
## |race : White|455|    68%  (156)   |    61%  (137)   |    64%  (293)   |
## +------------+---+-----------------+-----------------+-----------------+
## |    Black   |   |    25%  ( 58)   |    32%  ( 72)   |    29%  (130)   |
## +------------+---+-----------------+-----------------+-----------------+
## |    Other   |   |     7%  ( 17)   |     7%  ( 15)   |     7%  ( 32)   |
## +------------+---+-----------------+-----------------+-----------------+
## |apache      |454|  10.0/14.5/19.0 |  10.0/14.0/21.0 |  10.0/14.0/20.0 |
## +------------+---+-----------------+-----------------+-----------------+
## |temp0       |455| 99.9/100.8/101.7| 99.3/100.6/101.5| 99.5/100.7/101.6|
## +------------+---+-----------------+-----------------+-----------------+
## |temp1       |420| 99.2/100.5/101.5| 98.6/ 99.5/100.5| 98.9/ 99.9/101.1|
## +------------+---+-----------------+-----------------+-----------------+

L’évolution de la température (en °C) mesurée toutes les deux heures est résumée dans les deux figures suivantes, sous forme de données individuelles et agrégées sous forme de moyennes par période pour chacun des groupes de traitement.

p <- ggplot(data = na.omit(sepsis.long), aes(x = variable, y = value)) +
       geom_line(aes(group = id), col = grey(.7), alpha = .2) +
       geom_line(data = subset(na.omit(sepsis.long), variable %in% c("temp0", "temp1")),
                 aes(group = id, color = treat)) +
       scale_x_discrete(labels = seq(0, 6*2, by = 2)) +
       scale_y_continuous(breaks = seq(31.5, 43, by = 1.5)) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       labs(x = "Durée de suivi (heure)", y = "Température (°C)") +
       theme(legend.position = c(.5, 1.05), legend.direction = "horizontal")
p

d <- with(sepsis.long, summarize(value, llist(treat, variable), smean.cl.normal))

p <- ggplot(data = d, aes(x = variable, y = value, shape = treat, color = treat)) +
       geom_point(size = 2) +
       geom_line(aes(group = treat)) +
       geom_errorbar(aes(ymin=Lower, ymax=Upper), width=.1) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       scale_x_discrete(labels = seq(0, 6*2, by = 2)) +
       guides(color = FALSE, shape = FALSE) +
       geom_dl(aes(label = treat), method = list("smart.grid", cex = .8)) +
       labs(x = "Durée de suivi (heure)", y = "Température (°C)")
p

Enfin, on peut représenter la probabilité de survie au cours du temps à l’aide d’une courbe de Kaplan-Meier : en l’absence de censure, celle-ci représente directement la proportion de décès.

## 0 = alive, 1 = dead
st <- with(sepsis, Surv(time = followup, event = as.numeric(fate)-1))
s <- survfit(st ~ treat, data = sepsis)
p <- autoplot(s, censor = FALSE) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       scale_fill_manual("", values = c("steelblue", "orange")) +
       guides(fill = FALSE) +
       theme(legend.position = c(.9, .9)) +
       labs(x = "Followup time (hr.)", y = "Probability survival")
p

Approche par modèle mixte

Notons que la comparaison temp0 versus temp1 ne nécessiterait qu’une analyse de covariance. Ici, on souhaite modéliser l’ensemble des trajectoires individuelles dans les deux groupes. On aura donc quatre variables statistiques : l’identifiant du patient, le temps, le groupe de traitement et la réponse, et il faudra travailler avec le data frame en “format long” (sepsis.long). Le package lme4 est requis pour ces analyses (par commodité d’affichage et rapidité de calcul) mais on pourrait utiliser le package nlme également (la syntaxe correspondante est d’ailleurs indiquée en commentaire). Pour simplifier l’analyse, on ne se préoccupe pas des données manquantes (on ne garde que les patients avec les 7 mesures) et on ne considèrera que deux modèles simples :

cc <- aggregate(value ~ treat + id, sepsis.long, length)
idx <- cc$id[cc$value == 7]  ## complete case IDs
d <- subset(sepsis.long, id %in% idx)
d$id <- droplevels(d$id)     ## n = 328
names(d)[3] <- "time"
levels(d$time) <- 0:6
d$time <- as.numeric(as.character(d$time))

## nlme::lme(value ~ time * treat, data = d, random = ~ 1 | id)
m1 <- lmer(value ~ time * treat + (1 | id), data = d)
summary(m1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: value ~ time * treat + (1 | id)
##    Data: d
## 
## REML criterion at convergence: 5030
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.506 -0.508 -0.028  0.487  4.876 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.734    0.857   
##  Residual             0.351    0.592   
## Number of obs: 2296, groups:  id, 328
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)         37.99699    0.07264     523
## time                -0.07178    0.00858      -8
## treatIbuprofen      -0.23434    0.10467      -2
## time:treatIbuprofen -0.09926    0.01237      -8
## 
## Correlation of Fixed Effects:
##             (Intr) time   trtIbp
## time        -0.354              
## treatIbprfn -0.694  0.246       
## tm:trtIbprf  0.246 -0.694 -0.354

Le second modèle inclut une pente aléatoire pour chaque patient, en plus de l’intercept aléatoire :

## nlme::lme(value ~ time * treat, data = d, random = ~ time | id)
m2 <- lmer(value ~ time * treat + (time | id), data = d)

Les troisième et quatrième modèles incluent, respectivement, les deux mêmes effets aléatoires, mais cette fois sans corrélation, ou alors uniquement une pente aléatoire entre les deux groupes :

## nlme::lme(value ~ time * treat, data = d, random = list(id = pdDiag(~ time)))
m3 <- lmer(value ~ time * treat + (1 | id) + (0 + time | id), data = d)

## nlme::lme(value ~ time * treat, data = d, random = ~ 0 + time | id)
m4 <- lmer(value ~ time * treat + (0 + time | id), data = d)

screenreg(list(m1, m2, m3), include.variance = TRUE)

## 
## ==================================================================
##                           Model 1       Model 2       Model 3     
## ------------------------------------------------------------------
## (Intercept)                  38.00 ***     38.00 ***     38.00 ***
##                              (0.07)        (0.08)        (0.08)   
## time                         -0.07 ***     -0.07 ***     -0.07 ***
##                              (0.01)        (0.01)        (0.01)   
## treatIbuprofen               -0.23 *       -0.23 *       -0.23 *  
##                              (0.10)        (0.12)        (0.11)   
## time:treatIbuprofen          -0.10 ***     -0.10 ***     -0.10 ***
##                              (0.01)        (0.02)        (0.02)   
## ------------------------------------------------------------------
## AIC                        5042.04       4699.52       4767.09    
## BIC                        5076.47       4745.43       4807.26    
## Log Likelihood            -2515.02      -2341.76      -2376.54    
## Num. obs.                  2296          2296          2296       
## Num. groups: id             328           328           328       
## Var: id (Intercept)           0.73          1.01          0.87    
## Var: Residual                 0.35          0.23          0.24    
## Var: id time                                0.03                  
## Cov: id (Intercept) time                   -0.08                  
## Var: id.1 time                                            0.02    
## ==================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

Test du rapport de vraisemblance entre les trois premiers modèles

anova(m1, m3, m2)

## refitting model(s) with ML (instead of REML)

## Data: d
## Models:
## m1: value ~ time * treat + (1 | id)
## m3: value ~ time * treat + (1 | id) + (0 + time | id)
## m2: value ~ time * treat + (time | id)
##    Df  AIC  BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m1  6 5020 5054  -2504     5008                         
## m3  7 4747 4787  -2366     4733 274.85      1     <2e-16
## m2  8 4679 4725  -2332     4663  69.71      1     <2e-16

Approche par équations généralisées

Contrairement aux approches conditionnelles présentées plus haut, les équations généralisées permettent d’estimer l’évolution moyenne (et la structure de covariance) d’un ensemble d’observations corrélées. Les packages gee et geepack peuvent être utilisés, sachant que le package gee repose sur des erreurs-type asymptotiques. Les illustrations qui suivent reposent sur geepack.

Voici un modèle équivalent, sur le principe, au modèle à effets aléatoires précédents où l’on modélise l’évolution moyenne de la température en fonction du temps en interaction avec le groupe de traitement et en supposant une matrice de variance-covariance de travail dite “exchangeable” :

m <- geeglm(value ~ time * treat, data = d, id = id, corstr = "exch", scale.fix=TRUE)
summary(m)

## 
## Call:
## geeglm(formula = value ~ time * treat, data = d, id = id, corstr = "exch", 
##     scale.fix = TRUE)
## 
##  Coefficients:
##                     Estimate Std.err     Wald Pr(>|W|)
## (Intercept)          37.9970  0.0544 4.88e+05  < 2e-16
## time                 -0.0718  0.0153 2.21e+01  2.5e-06
## treatIbuprofen       -0.2343  0.0810 8.37e+00   0.0038
## time:treatIbuprofen  -0.0993  0.0221 2.01e+01  7.2e-06
## 
## Scale is fixed.
## 
## Correlation: Structure = exchangeable  Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha        0       0
## Number of clusters:   2296   Maximum cluster size: 1

Pour obtenir les matrices de variance-covariance estimée et de travail, il est nécessaire d’estimer les paramètres du modèle avec la fonction geese au lieu de geeglm :

m0 <- geese(value ~ time * treat, data = d, id = id, corstr = "exch", scale.fix=TRUE)
m0$vbeta

##           [,1]      [,2]      [,3]      [,4]
## [1,]  0.002959 -0.000685 -0.002959  0.000685
## [2,] -0.000685  0.000233  0.000685 -0.000233
## [3,] -0.002959  0.000685  0.006560 -0.001514
## [4,]  0.000685 -0.000233 -0.001514  0.000489

m0$vbeta.naiv

##           [,1]      [,2]      [,3]      [,4]
## [1,]  0.002949 -0.000681 -0.002949  0.000681
## [2,] -0.000681  0.000227  0.000681 -0.000227
## [3,] -0.002949  0.000681  0.006122 -0.001413
## [4,]  0.000681 -0.000227 -0.001413  0.000471

Enfin, voici les prédictions que donne cette approche marginale :

d$yhat3 <- predict(m)
p <- ggplot(data = d, aes(x = time, y = yhat3, color = treat)) +
       geom_line(aes(group = id), alpha = .2) +
       scale_x_continuous(breaks = 0:6, labels = seq(0, 6*2, by = 2)) +
       scale_y_continuous(breaks = seq(31.5, 43, by = 1.5)) +
       scale_color_manual("", values = c("steelblue", "orange")) +
       guides(color = FALSE) +
       geom_dl(aes(label = treat), method = "smart.grid") +
       labs(x = "Durée de suivi (heure)", y = "Température (°C)", caption = "Model GEE")
p

Exercice d’application

Les données proviennent d’une étude longitudinale de l’effet de la pollution atmosphérique sur l’état de santé mesuré par la capacité respiratoire des participants de l’étude. Les variables d’intérêt sont le statut respiratoire (resp) de 537 enfants âgés de 7 à 10 ans (age) et le statut fumeur de la mère (smoke).

data(ohio, package = "geepack")
str(ohio)

## 'data.frame':    2148 obs. of  4 variables:
##  $ resp : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ id   : int  0 0 0 0 1 1 1 1 2 2 ...
##  $ age  : int  -2 -1 0 1 -2 -1 0 1 -2 -1 ...
##  $ smoke: int  0 0 0 0 0 0 0 0 0 0 ...

ohio$age <- ohio$age + 9    ## in years

Le tableau suivant résume la proportion de cas en fonction de l’âge (lignes) et du statut fumeur (colonne) :

summary(resp ~ age + smoke, data = ohio, method = "cross")

## 
##  mean by age, smoke 
## 
## +----+
## |N   |
## |resp|
## +----+
## +---+-----+-----+-----+
## |age|  0  |  1  | ALL |
## +---+-----+-----+-----+
## |7  | 350 | 187 | 537 |
## |   |0.160|0.166|0.162|
## +---+-----+-----+-----+
## |8  | 350 | 187 | 537 |
## |   |0.149|0.209|0.169|
## +---+-----+-----+-----+
## |9  | 350 | 187 | 537 |
## |   |0.143|0.187|0.158|
## +---+-----+-----+-----+
## |10 | 350 | 187 | 537 |
## |   |0.106|0.139|0.117|
## +---+-----+-----+-----+
## |ALL|1400 | 748 |2148 |
## |   |0.139|0.175|0.152|
## +---+-----+-----+-----+

Estimer les paramètres d’un modèle GEE avec une matrice de variance-covariance “exchangeable” et un terme d’interaction age:smoke. Attention, la variable étant binaire, il est nécessaire de préciser une famille de type binomial. Le statut fumeur de la mère est-il significatif au seuil conventionel de 5 % ? Même question pour l’âge.
Calculer l’odds-ratio ajusté pour la variable age, avec son son intervalle de confiance à 95 %.
Représenter graphiquement les prédictions moyennes de ce modèle.
Comparer avec une approche par modèle mixte.

Méthodes biostatistiques

Modèles pour données corrélées

Modèles à effets aléatoires

Chargement des données

Description des variables

Modélisation

Cas des données longitudinales

Chargement des données

Description des variables

Approche par modèle mixte

Approche par équations généralisées

Exercice d’application