Environnement et packages

# activer l'environnement
using Pkg
Pkg.activate("env_julia_regression")

  Activating project at `c:\Users\ricco\Desktop\demo\env_julia_regression`

# liste des packages installés
Pkg.status()

Status `C:\Users\ricco\Desktop\demo\env_julia_regression\Project.toml`
  [336ed68f] CSV v0.10.16
  [13f3f980] CairoMakie v0.15.11
  [a93c6f00] DataFrames v1.8.2
  [31c24e10] Distributions v0.25.126
  [38e38edf] GLM v1.9.5
  [7073ff75] IJulia v1.34.4
  [43a3c2be] PairPlots v3.0.3
  [f3b207a7] StatsPlots v0.15.8

Importation des données

# packages
import DataFrames as DFR
import CSV

# lecture des données
df_all = CSV.read("vehicules_1.txt", DFR.DataFrame ; delim = "\t")

# premières lignes
DFR.first(df_all,5)

5×5 DataFrame

Row	modele	cylindree	puissance	poids	conso
	String31	Int64	Int64	Int64	Float64
1	Maserati Ghibli GT	2789	209	1485	14.5
2	Daihatsu Cuore	846	32	650	5.7
3	Toyota Corolla	1331	55	1010	7.1
4	Fort Escort 1.4i PT	1390	54	1110	8.6
5	Mazda Hachtback V	2497	122	1330	10.8

# isoler les variables actives
df = DFR.select(df_all,DFR.Not(:modele))

# description
# on obtient notamment le type interne des variables (eltype)
DFR.describe(df)

4×7 DataFrame

Row	variable	mean	min	median	max	nmissing	eltype
	Symbol	Float64	Real	Float64	Real	Int64	DataType
1	cylindree	2028.64	846	1984.0	5987	0	Int64
2	puissance	91.44	29	85.0	300	0	Int64
3	poids	1257.8	650	1155.0	2250	0	Int64
4	conso	9.648	5.7	9.0	18.7	0	Float64

# dimensions
println(DFR.size(df))

# nombre d'observations
n = DFR.nrow(df)
println(n)

# nombre de variables explicatives
p = DFR.ncol(df) - 1
println(p)

(25, 4)
25
3

# inspection rapide des données - boxplot avec StatsPlots
import StatsPlots as STP
# un boxplot par variable
p1 = STP.boxplot(df.cylindree,title="Cylindree",label="")
p2 = STP.boxplot(df.puissance,title="Puissance",label="")
p3 = STP.boxplot(df.poids,title="Poids",label="")
p4 = STP.boxplot(df.conso,title="Conso",label="")

# affichage dans une grille (2 x 2)
STP.plot(p1,p2,p3,p4,layout=(2,2),size=(600,600))

# et un pairplot (PairPlots)
import PairPlots as PP
import CairoMakie as CM
PP.pairplot(df => (PP.Scatter(markersize=8),
                    PP.MarginHist(),
                    PP.PearsonCorrelation()
                   ))

Régression linéaire multiple

Lancement des calculs

# importer la librairie
import GLM

# définir et lancer la régression
modele = GLM.lm(GLM.@formula(conso ~ cylindree + puissance + poids), df)

# affichage
display(modele)

StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

conso ~ 1 + cylindree + puissance + poids

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.   Std. Error      t  Pr(>|t|)    Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)   1.29984      0.614073      2.12    0.0464   0.0228017   2.57687
cylindree    -0.000473994  0.000497236  -0.95    0.3513  -0.00150805  0.000560066
puissance     0.0302691    0.0074708     4.05    0.0006   0.0147327   0.0458054
poids         0.00520108   0.000806704   6.45    <1e-05   0.00352345  0.00687872
─────────────────────────────────────────────────────────────────────────────────

Accès aux résultats détaillés

# informations

# R2
println("R2 = $(GLM.r2(modele))")
# R2 ajusté
println("R2-adj = $(GLM.adjr2(modele))")
# AIC
println("AIC = $(GLM.aic(modele))")
# BIC
println("BIC = $(GLM.bic(modele))")
# degrés de liberté des résidus
println("ddl modele = $(GLM.dof_residual(modele))")

R2 = 0.9567033306282597
R2-adj = 0.9505180921465826
AIC = 58.80056270070373
BIC = 64.89494182504473
ddl modele = 21.0

# SCT
println("SCT = $(GLM.nulldeviance(modele))")
# SCR
println("SCR = $(GLM.deviance(modele))")
# Test F de significativité globale
dev_modele = GLM.deviance(modele)
dev_nulle = GLM.nulldeviance(modele)
FTest = ((dev_nulle - dev_modele)/p)/(dev_modele/GLM.dof_residual(modele))
println("F satistic = $FTest")

SCT = 238.10239999999996
SCR = 10.30904088941784
F satistic = 154.6752535835679

# et pour obtenir la p-value
import Distributions as distrib
p_value = distrib.ccdf(distrib.FDist(p,GLM.dof_residual(modele)),FTest)
println("p_value = $p_value")

p_value = 1.7868899968536525e-14

Comparaison de modèles

# régression avec puissance et poids
modele_simple = GLM.lm(GLM.@formula(conso ~ puissance + poids), df)
display(modele_simple)

StatsModels.TableRegressionModel{GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

conso ~ 1 + puissance + poids

Coefficients:
────────────────────────────────────────────────────────────────────────────
                  Coef.   Std. Error     t  Pr(>|t|)   Lower 95%   Upper 95%
────────────────────────────────────────────────────────────────────────────
(Intercept)  1.3593      0.609627     2.23    0.0363  0.0950152   2.62359
puissance    0.0244309   0.00426939   5.72    <1e-05  0.0155767   0.033285
poids        0.00481375  0.000695453  6.92    <1e-06  0.00337147  0.00625603
────────────────────────────────────────────────────────────────────────────

# f-test de comparaison des modèles imbriqués
# nullité du coefficient de cylindree
GLM.ftest(modele_simple.model,modele.model)

F-test: 2 models fitted on 25 observations
────────────────────────────────────────────────────────────────
     DOF  ΔDOF      SSR     ΔSSR      R²     ΔR²      F*   p(>F)
────────────────────────────────────────────────────────────────
[1]    4        10.7551           0.9548                        
[2]    5     1  10.3090  -0.4461  0.9567  0.0019  0.9087  0.3513
────────────────────────────────────────────────────────────────

Diagnostic de la régression

Graphique des résidus

# accès aux résidus
residus = GLM.residuals(modele)
# somme des carrés = SCR
sum(residus .^ 2)

10.30904088941784

# graphique des résidus (CairoMakie)
fig = CM.Figure()
ax = CM.Axis(fig[1,1],title = "Graphique des résidus")
CM.scatter!(df.conso,residus,markersize=7,color=:gray)
CM.hlines!([0],linewidth=2,color=:blue)
display(fig)

CairoMakie.Screen{IMAGE}

Points atypiques et influents

# matrice des explicatives incluant la constante
X = GLM.modelmatrix(modele)
X

25×4 Matrix{Float64}:
 1.0  2789.0  209.0  1485.0
 1.0   846.0   32.0   650.0
 1.0  1331.0   55.0  1010.0
 1.0  1390.0   54.0  1110.0
 1.0  2497.0  122.0  1330.0
 1.0  2473.0  125.0  1570.0
 1.0  2438.0   97.0  1800.0
 1.0  2165.0  101.0  1500.0
 1.0  1396.0   66.0  1140.0
 1.0  1984.0   85.0  1155.0
 ⋮                   
 1.0  1195.0   33.0   895.0
 1.0  1597.0   74.0  1080.0
 1.0  1761.0   74.0  1100.0
 1.0  1998.0   66.0  1300.0
 1.0  1998.0   89.0  1140.0
 1.0  1998.0   89.0  1560.0
 1.0   899.0   29.0   730.0
 1.0  1984.0   85.0  1635.0
 1.0  1242.0   55.0   940.0

# calcul (X'X)(-1)
import LinearAlgebra as LA
XtX_inv = LA.inv(X' * X)
XtX_inv

4×4 Matrix{Float64}:
  0.768142      6.31883e-5   0.00213739   -0.000836199
  6.31883e-5    5.03648e-7  -6.20345e-6   -4.11564e-7
  0.00213739   -6.20345e-6   0.000113694   4.05781e-8
 -0.000836199  -4.11564e-7   4.05781e-8    1.32565e-6

Levier

# calcul du levier [diagonale de Hat matrix]
# (pas de fonction directe dans GLM)
h = [X[i,:]' * XtX_inv * X[i,:] for i in axes(X,1)]
h

25-element Vector{Float64}:
 0.7218393114156288
 0.17493952912740662
 0.060524388382147576
 0.059843026531024676
 0.058334574215776674
 0.09834564668360823
 0.30693840415549445
 0.09434880660512285
 0.07277268588656306
 0.052437822052703476
 ⋮
 0.10112178818413445
 0.05400760473146257
 0.051215731641375425
 0.10772156064633785
 0.05567033097157892
 0.16884878723671604
 0.1320403074389517
 0.2444283179546619
 0.07603256793299977

# point suspect si levier > (2(p+1))/n
seuil_h = 2.0*(p+1)/n
seuil_h

0.32

# véhicules avec leviers élevés
df_all.modele[h .> seuil_h]

2-element Vector{InlineStrings.String31}:
 "Maserati Ghibli GT"
 "Mercedes S 600"

Résidus standardisés

# écart-type des résidus la régression
# sqrt() de SCR / degrés de libertés
sigma_err = sqrt(GLM.deviance(modele)/GLM.dof_residual(modele))
sigma_err

0.7006473499699383

# résidus standardisés
r_stand = residus ./ (sigma_err .* sqrt.(1.0 .- h))
r_stand

25-element Vector{Float64}:
  1.2780867370612505
  0.7099857724709607
 -0.7168816533538167
  0.8114738324543007
  0.10804304213931906
  0.9364392287014734
  0.6132625750534932
  0.851199122920511
 -1.282304557644014
  0.8217310425248503
  ⋮
  0.6213315630020293
 -1.4673536176653375
  0.8407170734490949
 -2.2845031425763938
 -0.2584637430707878
 -0.5642691516935364
  0.8451731756954376
  0.2691584061811583
 -0.9873325108997584

Résidus studentisés

# résidus studentisés à partir
# des résidus standardisés et des leviers
r_stud = r_stand .* sqrt.((n-p-2)./((n-p-1).-r_stand.^2))
r_stud

25-element Vector{Float64}:
  1.2988225490017054
  0.7013437529280638
 -0.7083257403648012
  0.804633149800432
  0.10546853024876957
  0.9335710001268489
  0.6039152104492624
  0.845397203104742
 -1.3034722759858302
  0.8151396128568138
  ⋮
  0.6120090372666391
 -1.5115764906124813
  0.8346216631144779
 -2.5718099554542344
 -0.2526369399270822
 -0.5548929864068488
  0.8392008070906738
  0.26312596913356173
 -0.9867117065288005

# calculer le seuil avec la loi de Student
# quantile à 95% bilatéral (c.-à-d. à 0.975)
seuil_r_stud = distrib.quantile(distrib.TDist(n-p-2),0.975)
seuil_r_stud

2.0859634472658644

# atypiques au sens des résidus studentisé
df_all.modele[abs.(r_stud) .> seuil_r_stud]

2-element Vector{InlineStrings.String31}:
 "Hyundai Sonata 3000"
 "Mitsubishi Galant"

# visualiser dans le graphique des résidus
fig = CM.Figure()
ax = CM.Axis(fig[1,1],title = "Graphique des résidus studentisés")
CM.scatter!(df.conso,r_stud,markersize=7,color=:gray)
CM.hlines!([0],linewidth=1,color=:blue)
CM.hlines!([-seuil_r_stud],linewidth=2,color=:red,linestyle=:dash)
CM.hlines!([+seuil_r_stud],linewidth=2,color=:red,linestyle=:dash)
display(fig)

CairoMakie.Screen{IMAGE}

# même graphique en faisant apparaître le nom des modèles
fig = CM.Figure()
ax = CM.Axis(fig[1,1],title = "Graphique des résidus studentisés")
CM.scatter!(df.conso,r_stud,markersize=7,color=:silver)
CM.hlines!([0],linewidth=1,color=:blue)
CM.hlines!([-seuil_r_stud],linewidth=2,color=:orange,linestyle=:dash)
CM.hlines!([+seuil_r_stud],linewidth=2,color=:orange,linestyle=:dash)
# pour l'ensemble des voitures suspectes
for voiture in df_all.modele[abs.(r_stud) .> seuil_r_stud]
    id = findfirst(==(voiture),df_all.modele)
    CM.text!([df.conso[id]],[r_stud[id]];
            text=voiture,color=:red,align=(:center,:center))    
end
display(fig)

CairoMakie.Screen{IMAGE}

Points observés vs. prédits

# prédiction en resubstitution
# on ne passe pas les données à traiter
# la fonction utilise automatiquement les données d'apprentissage
pred_resub = GLM.predict(modele)
pred_resub

25-element Vector{Float64}:
 14.027711154840901
  5.24815235845113
  7.586843868731834
  8.048717521678618
 10.726541069633335
 12.076984108874377
 12.442289356044935
 11.13244084551341
  8.565134792067665
  8.939554952087917
  ⋮
  6.387263072946842
  8.399949559541943
  8.42623627353742
  9.111963973351832
  8.975979034452104
 11.160434051370197
  5.548310183822661
 11.43607497142288
  7.264953464648552

# graphique y observé vs. y predit
# min et max pour un graphique "carré"
min_val = minimum([minimum(df.conso),minimum(pred_resub)])
max_val = maximum([maximum(df.conso),maximum(pred_resub)])

# graphique
fig = CM.Figure()
ax = CM.Axis(fig[1,1],xlabel="Valeurs Observées",ylabel="Valeurs Prédites")
CM.xlims!(ax,min_val,max_val+1)
CM.ylims!(ax,min_val,max_val+1)
CM.scatter!(df.conso,pred_resub)
CM.ablines!(0,1,color=:green)
display(fig)

CairoMakie.Screen{IMAGE}

Prédiction ponctuelle et par intervalle

Importation des observations supplémentaires

# lecture des données
df_supp = CSV.read("vehicules_2.txt", DFR.DataFrame ; delim = "\t")
df_supp

6×5 DataFrame

Row	modele	cylindree	puissance	poids	conso
	String31	Int64	Int64	Int64	Float64
1	Nissan Primera 2.0	1997	92	1240	9.2
2	Fiat Tempra 1.6 Liberty	1580	65	1080	9.3
3	Opel Omega 2.5i V6	2496	125	1670	11.3
4	Subaru Vivio 4WD	658	32	740	6.8
5	Seat Ibiza 2.0 GTI	1983	85	1075	9.5
6	Ferrari 456 GT	5474	325	1690	21.3

# variables explicatives
X_supp = df_supp[:,[:cylindree,:puissance,:poids]]
X_supp

6×3 DataFrame

Row	cylindree	puissance	poids
	Int64	Int64	Int64
1	1997	92	1240
2	1580	65	1080
3	2496	125	1670
4	658	32	740
5	1983	85	1075
6	5474	325	1690

# nombre d'observation
n_supp = size(X_supp,1)
n_supp

6

Prédiction ponctuelle

# prédiction ponctuelle
# on passe cette fois-ci la matrice des descripteurs
pred_supp = GLM.predict(modele,X_supp)
pred_supp

6-element Vector{Union{Missing, Float64}}:
  9.58736854815561
  8.135585902890996
 12.586190593004323
  5.805360662403514
  8.523942275817518
 17.332471445359175

# confrontation points observés vs. prédits
# min et max pour un graphique "carré"
min_val = minimum([minimum(df_supp.conso),minimum(pred_supp)])
max_val = maximum([maximum(df_supp.conso),maximum(pred_supp)])

# graphique
fig = CM.Figure()
ax = CM.Axis(fig[1,1],xlabel="Valeurs Observées",ylabel="Valeurs Prédites")
CM.xlims!(ax,min_val-1,max_val+1)
CM.ylims!(ax,min_val-1,max_val+1)
CM.scatter!(df_supp.conso,pred_supp,color=:blue)
CM.ablines!(0,1,color=:green)
display(fig)

CairoMakie.Screen{IMAGE}

Variances de l'erreur de prédiction

# ajouter une constante 1 à la matrice des X_supp
X1_supp = hcat(ones(1,n_supp)',Matrix(X_supp))
X1_supp

6×4 Matrix{Float64}:
 1.0  1997.0   92.0  1240.0
 1.0  1580.0   65.0  1080.0
 1.0  2496.0  125.0  1670.0
 1.0   658.0   32.0   740.0
 1.0  1983.0   85.0  1075.0
 1.0  5474.0  325.0  1690.0

# rappel - écart-type de la régression
sigma_err

0.7006473499699383

# rappel (X'X)(-1)
XtX_inv

4×4 Matrix{Float64}:
  0.768142      6.31883e-5   0.00213739   -0.000836199
  6.31883e-5    5.03648e-7  -6.20345e-6   -4.11564e-7
  0.00213739   -6.20345e-6   0.000113694   4.05781e-8
 -0.000836199  -4.11564e-7   4.05781e-8    1.32565e-6

# variance de l'erreur de prédiction pour chaque individu
var_err = zeros(n_supp)

# pour chaque individu à traiter
for i in 1:n_supp
    # description - vecteur à caster en matrice
    tmp = reshape(X1_supp[i,:],1,:)
    # produit matriciel
    pm = tmp * XtX_inv * tmp'
    # matrice (1,1) -> transformer en scalaire
    pm = pm[1,1]
    # variance
    var_err[i] = sigma_err^2 * (1.0 + pm)
end

# affichage
var_err

6-element Vector{Float64}:
 0.5108941292423475
 0.5156051122226518
 0.5651575907827666
 0.5649406208483717
 0.5300043071227483
 1.1128200339955032

Intervalle de confiance de prédiction

# quantile de la loi de Student
# pour un niveau de confiance à 95%
qt = distrib.quantile(distrib.TDist(n-p-1),0.975)
qt

2.07961384472768

# bornes basses
yb = pred_supp .- qt .* sqrt.(var_err)
yb

6-element Vector{Float64}:
  8.100925904265084
  6.642305700352399
 11.02280004278032
  4.242270242013515
  7.009954392344032
 15.138680862757342

# bornes hautes
yh = pred_supp .+ qt .* sqrt.(var_err)
yh

6-element Vector{Float64}:
 11.073811192046135
  9.628866105429594
 14.149581143228326
  7.368451082793512
 10.037930159291005
 19.52626202796101

# graphique - intervalle de confiance de prédiction
# vs. valeurs observées de la cible
# pour les individus supplémentaires
fig = CM.Figure()
ax = CM.Axis(fig[1,1],
            xticks=(1:n_supp,df_supp.modele),
            xticklabelrotation=π/2)
CM.rangebars!(ax,1:n_supp,yb,yh,whiskerwidth=10)
CM.scatter!(ax,1:n_supp,df_supp.conso,color=:red)
display(fig)

CairoMakie.Screen{IMAGE}