Diagnosticarea modelului: statistici vs învățare automată

În această postare, arătăm cât de diferit cazuri de utilizare necesită diferit Diagnostic de model. Pe scurt, comparăm (statistic) Inferență şi predicție.

Ca exemplu, folosim un model liniar simplu pentru setul de date al indicelui Munchen Rent, care a fost furnizat cu drag de autorii regresiei – modele, metode și aplicații a doua ed. (2021). Acest set de date conține chirii lungi în EUR (rent) pentru aproximativ 3000 de apartamente la Munchen, Germania, din 1999. Apartamentele au mai multe caracteristici precum zona de zi (area) în metri pătrați, anul construcției (yearc), calitatea locației (location0: medie, 1: bun, 2: sus), calitatea camerelor de baie (bath0: standard, 1: premium), calitatea bucătăriei (kitchen0: standard, 1: premium), indicator pentru încălzirea centrală (cheating)

Variabila țintă este Y=text{rent} Iar obiectivul modelului nostru este de a prezice chiria medie, E(Y) (omitem condiționarea pe X pentru brevet).

Disclaimer: Înainte de a prezenta cazurile de utilizare, permiteți -mi să afirm clar că nu sunt în afacerea cu închirieri de apartamente și totul aici este doar în scopul de a demonstra o bună practică statistică.

Inferență

Primul caz de utilizare este despre inferența efectului caracteristicilor. Imaginează -ți punctul de vedere al unui investitor care dorește să știe dacă instalarea unei încălziri centrale merită (financiar). Pentru a pune terenul pe care să se bazeze o decizie, un statisticist trebuie să aibă răspunsuri la:

• Ce este efect a variabilei cheating pe chirie.
• Este acest efect semnificativ statistic?

Predicție

Al doilea caz de utilizare este despre predicție. De această dată, luăm punctul de vedere al cuiva care caută un apartament nou de închiriat. Pentru a ști dacă chiria propusă de către proprietar este aproximativ corectă sau necorespunzătoare (prea mare), o valoare de referință ar fi foarte convenabilă. Se poate cere vecinilor sau se poate cere unui model să prezică chiria apartamentului în cauză.

Model potrivire

Înainte de a răspunde la întrebările de mai sus și de a face unele diagnostice cheie, trebuie să încărcăm datele și să ne potrivim unui model. Alegem un model liniar simplu și un model direct rent.

Note:

• Pentru indicii de închiriere, precum și pentru prețurile locuințelor, de multe ori se transformă log-ul variabilului țintă înainte de modelare sau unul folosește o legătură de jurnal și o funcție de pierdere adecvată (de exemplu, devianța gamma).
• Versiunea noastră Python folosește GeneralizedLinearRegressor din pachet glum. La fel de bine am fi ales alte implementări precum StatsModels.Regression.Linear_Model.Ols. În acest fel, trebuie să implementăm singuri diagnosticul rezidual, ceea ce face clar ce este reprezentat.

Pentru brevet, omitem importurile și încărcarea datelor. Modelul nostru este apoi potrivit de:

lm = glum.GeneralizedLinearRegressor(
    alpha=0,
    drop_first=True,  # this is very important if alpha=0
    formula="bs(area, degree=3, df=4) + yearc"
      	" + C(location) + C(bath) + C(kitchen) + C(cheating)"
)
lm.fit(X_train, y_train)

model = lm(
  formula = rent ~ bs(area, degree = 3, df = 4) + yearc + location + bath + kitchen + cheating,
  data = df_train
)

Diagnostice pentru inferență

Tabelul de coeficient ne va spune deja efectul cheating variabil. Pentru modele mai implicate, cum ar fi copacii impulsionați de gradient sau plase neuronale, se pot folosi dependență parțială și valori de formă pentru a evalua efectul caracteristicilor.

lm.coef_table(X_train, y_train)

summary(model)
confint(model)

Vedem asta Ceteris paribusceea ce înseamnă că toate celelalte egale, o încălzire centrală crește chiria lunară cu aproximativ 108 EUR. Nu dimensiunea efectului de 108 EUR, dar faptul că există un efect al încălzirii centrale asupra chirii pare semnificativ statistic:
Acest lucru este indicat de probabilitatea foarte mică, adică valoarea p, pentru hipoteza nulă a cheating având un coeficient de zero.
De asemenea, vedem că intervalul de încredere cu nivelul implicit de încredere de 95%: (ci_lower, ci_upper) = (87, 129).
Acest lucru arată incertitudinea efectului estimat.

Pentru o clădire cu 10 apartamente și cu un orizont de investiții de aproximativ 10 ani, efectul estimat oferă aproximativ un buget de 13000 EUR (intervalul este de aproximativ 10500 până la 15500 cu încredere de 95%).

Un statisticist bun ar trebui să pună mai multe întrebări suplimentare:

• Setul de date este la îndemână o bună reprezentare a populației?
• Există confuzii sau efecte de interacțiune, în special între cheating Și alte caracteristici?
• Sunt ipotezele pentru valoarea p scăzută și intervalul de încredere al cheating valabil?

Aici, vom aborda doar ultima întrebare și chiar și cea parțială. Ce presupuneri au fost făcute? Termenul de eroare, epsilon = Y - E(Y)ar trebui să fie homoscedastic și distribuit normal. Deoarece eroarea nu este observabilă (deoarece adevărat model pentru E(Y) este necunoscut), unul înlocuiește E(Y) prin predicția modelului hat{E}(Y)acest lucru dă reziduurilor, hat{epsilon} = Y - hat{E}(Y) = y - text{fitted values}în schimb. Pentru homoscedasticitate, reziduurile ar trebui să arate ca zgomot alb (aleatoriu). Normalitatea, pe de altă parte, devine mai puțin o preocupare cu datele mai mari datorită teoremei limită centrală. Cu aproximativ 3000 de puncte de date, suntem departe de date micidar ar putea fi totuși o idee bună să verificați normalitatea.

Instrumentele de diagnostic pentru a verifica dacă sunt parcele reziduale și cuantile-quatile (QQ).

# See notebook for a definition of residual_plot.
import seaborn as sns
fig, axes = plt.subplots(ncols=2, figsize=(4.8 * 2.1, 6.4))
ax = residual_plot(model=lm, X=X_train, y=y_train, ax=axes(0))
sns.kdeplot(
    x=lm.predict(X_train),
    y=residuals(lm, X_train, y_train, kind="studentized"),
    thresh=.02,
    fill=True,
    ax=axes(1),
).set(
    xlabel="fitted",
    ylabel="studentized residuals",
    title="Contour Plot of Residuals",
)

autoplot(model, which = c(1, 2))  # from library ggfortify
# density plot of residuals
ggplot(model, aes(x = .fitted, y = .resid)) + geom_point() +
  geom_density_2d() + geom_density_2d_filled(alpha = 0.5)

Cu cât este mai multe puncte de date, cu atât mai puțin informativ este un complot de împrăștiere. Prin urmare, punem un complot de contur pe dreapta.

Perspective vizuale:

• Se pare că există o variabilitate mai mare pentru valori mai mari montate. Acesta este un indiciu că homoscedasticitatea ar putea fi încălcată.
• Reziduurile par a fi centrate în jurul orei 0. Acesta este un indiciu că modelul este bine calibrat (adecvat).

# See notebook for a definition of qq_plot.
qq_plot(lm, X_train, y_train)

autoplot(model, which = 2)

Plotul QQ arată cuantile distribuției teoretice asumate a reziduurilor pe axa X și valorile ordonate ale reziduurilor pe axa y. În versiunea Python, am decis să folosim reziduurile de student, deoarece normalitatea erorii implică o distribuție a studentului (T) pentru aceste reziduuri.

Observații finale:

• Am putea face parcele similare pe eșantionul de testare, dar nu avem nevoie neapărat de un eșantion de testare pentru a răspunde la întrebările de inferență.
• Este o practică bună să complotați reziduurile față de fiecare dintre caracteristici.

Diagnostic pentru predicție

Dacă suntem interesați doar de predicțiile chirii medii, hat{E}(Y)nu ne pasă prea mult de distribuția probabilității Y. Vrem doar să știm dacă predicțiile sunt suficient de aproape de media reală a chirii E(Y). Într -un argument similar ca pentru termenul de eroare și reziduuri, trebuie să acceptăm acest lucru E(Y) nu este observabil (este cantitatea pe care vrem să o prezicem). Deci trebuie să ne întoarcem la observațiile Y Pentru a judeca dacă modelul nostru este bine calibrat, adică închideți idealul E(Y).

Foarte important, aici folosim eșantionul de testare în toate diagnosticele noastre, deoarece Ne temem de prejudecata din eșantion.

Începem simplu printr -o privire asupra calibrării necondiționate, adică reziduul mediu (negativ) frac{1}{n}sum(hat{E}(Y_i)-Y_i).

compute_bias(
    y_obs=np.concatenate((y_train, y_test)),
    y_pred=lm.predict(pd.concat((X_train, X_test))),
    feature=np.array(("train") * X_train.shape(0) + ("test") * X_test.shape(0)),
)

print(paste("Train set mean residual:", mean(resid(model))))
print(paste("Test set mean residual: ", mean(df_test$rent - predict(model, df_test))))

Nu este de mirare că bias_mean În setul de tren este aproape zero.
Acesta este Proprietate de echilibru de modele liniare (generalizate) (cu termen de interceptare). Cu toate acestea, pe setul de testare, detectăm în medie o mică prejudecată de aproximativ 2 EUR pe apartament.

În continuare, avem un aspect diagrame de fiabilitate care conțin mult mai multe informații despre calibrarea și prejudecata unui model decât calibrarea necondiționată de mai sus. De fapt, evaluează auto-calibrarea, adică cât de bine folosește modelul propriilor informații.
Un model ideal s -ar afla pe linia diagonală punctată.

fig, axes = plt.subplots(ncols=2, figsize=(4.8 * 2.1, 6.4))
plot_reliability_diagram(y_obs=y_train, y_pred=lm.predict(X_train), n_bootstrap=100, ax=axes(0))
axes(0).set_title(axes(0).get_title() + f" train set (n={X_train.shape(0)})")
plot_reliability_diagram(y_obs=y_test, y_pred=lm.predict(X_test), n_bootstrap=100, ax=axes(1))
axes(1).set_title(axes(1).get_title() + f" test set (n={X_test.shape(0)})")

iso_train = isoreg(x = model$fitted.values, y = df_train$rent)
iso_test = isoreg(x = predict(model, df_test), y = df_test$rent)
bind_rows(
  tibble(set = "train", x = iso_train$x(iso_train$ord), y = iso_train$yf),
  tibble(set = "test", x = iso_test$x(iso_test$ord), y = iso_test$yf),
) |>
  ggplot(aes(x=x, y=y, color=set)) + geom_line() +
  geom_abline(intercept = 0, slope = 1, linetype="dashed") +
  ggtitle("Reliability Diagram")

Perspective vizuale:

• Graficele de pe tren și set de testare arată foarte asemănător.
  Intervalele de incertitudine mai mari de pe setul de testare provin din faptul că este o dimensiune mai mică a eșantionului.
• Modelul pare să se întindă în jurul diagonalei care indică o auto-calibrare bună pentru cea mai mare parte a gamei.
• Valorile prezise foarte mari par a fi în mod sistematic prea scăzute, adică graficul este deasupra diagonalei.

În cele din urmă, evaluăm calibrarea condiționată, adică calibrarea cu privire la caracteristici. Prin urmare, complotăm unul dintre graficele noastre preferate pentru fiecare caracteristică. Este format din:

• Valoarea medie observată a Y pentru fiecare valoare (bazată) a caracteristicii
• Valoarea medie prevăzută
• dependență parțială
• Histograma caracteristicii (axa y dreaptă)

fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(12, 5*4), sharey=True)
for i, col in enumerate(("area", "yearc", "bath", "kitchen", "cheating")):
    plot_marginal(
        y_obs=y_train,
        y_pred=lm.predict(X_train),
        X=X_train,
        feature_name=col,
        predict_function=lm.predict,
        ax=axes(i)(0),
    )
    plot_marginal(
        y_obs=y_test,
        y_pred=lm.predict(X_test),
        X=X_test,
        feature_name=col,
        predict_function=lm.predict,
        ax=axes(i)(1),
    )
    axes(i)(0).set_title("Train")
    axes(i)(1).set_title("Test")
    if i != 0:
        axes(i)(0).get_legend().remove()
    axes(i)(1).get_legend().remove()
fig.tight_layout()

xvars = c("area", "yearc", "bath", "kitchen", "cheating")
m_train = feature_effects(model, v = xvars, data = df_train, y = df_train$rent)
m_test = feature_effects(model, v = xvars, data = df_test, y = df_test$rent)

c(m_train, m_test) |> 
  plot(
    share_y = "rows",
    ncol = 2,
    byrow = FALSE,
    stats = c("y_mean", "pred_mean", "pd"),
    subplot_titles = FALSE,
    # plotly = TRUE,
    title = "Left: Train - Right: Test",
  )

Perspective vizuale:

• Pe setul de trenuri, caracteristicile categorice par să aibă o calibrare perfectă, deoarece media observată este egală în medie. Acesta este din nou un rezultat al proprietății soldului. În setul de testare, vedem o abatere, în special pentru nivelul categoric cu o dimensiune mai mică a eșantionului. Aceasta este o demonstrație bună de ce complotarea atât a trenului, cât și a setului de testare este o idee bună.
• Zona de caracteristici numerice și anul construcției par bine, dar o privire mai atentă nu poate face rău.

Următorul efectuăm un complot de prejudecăți, care complotează diferența medie de minus prevăzută observată pe valoarea caracteristicii. Valorile ar trebui să fie în jur de zero, astfel încât să putem face zoom pe axa y.
Acest lucru este foarte asemănător cu complotul rezidual, dar informațiile sunt mai bine condensate pentru scopul său.

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 2*4), sharey=True)
axes(0,0).set_ylim(-150, 150)
for i, col in enumerate(("area", "yearc")):
    plot_bias(
        y_obs=y_train,
        y_pred=lm.predict(X_train),
        feature=X_train(col),
        ax=axes(i)(0),
    )
    plot_bias(
        y_obs=y_test,
        y_pred=lm.predict(X_test),
        feature=X_test(col),
        ax=axes(i)(1),
    )
    axes(i)(0).set_title("Train")
    axes(i)(1).set_title("Test")
fig.tight_layout()

c(m_train(c("area", "yearc")), m_test(c("area", "yearc"))) |> 
  plot(
    ylim = c(-150, 150),
    ncol = 2,
    byrow = FALSE,
    stats = "resid_mean",
    subplot_titles = FALSE,
    title = "Left: Train - Right: Test",
    # plotly = TRUE,
    interval = "ci"
  )

Perspective vizuale:

• Pentru valori mari de area şi yearc În anii ’40 -’50, există doar puține observații disponibile. Totuși, modelul ar putea fi îmbunătățit pentru acele regiuni.
• Prejudecata yearc prezintă o curbă parabolică. Efectul liniar simplu în modelul nostru pare prea simplist. Un model rafinat ar putea folosi splines în schimb, în ​​ceea ce privește area.

Observații finale:

• Predicțiile pentru area Mai mare decât în ​​jur de 120 de metri pătrați și pentru anul de construcție în jurul celui de -al doilea război mondial sunt mai puțin de încredere.
• Pentru toate celelalte, prejudecata este mai mică de 50 EUR în medie.
  Prin urmare, aceasta este o estimare brută a incertitudinii de predicție.
  Ar trebui să fie suficient pentru a preveni chirii necorespunzătoare (sau mici) (în medie).

Codul complet Python și R este disponibil sub: