[Appendice — Regressione multipla

# Standard library imports
import os

# Third-party imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import arviz as az
import scipy.stats as stats
from scipy.special import expit  # Funzione logistica
from cmdstanpy import cmdstan_path, CmdStanModel

# Configuration
seed = sum(map(ord, "stan_poisson_regression"))
rng = np.random.default_rng(seed=seed)
az.style.use("arviz-darkgrid")
%config InlineBackend.figure_format = "retina"

# Define directories
home_directory = os.path.expanduser("~")
project_directory = f"{home_directory}/_repositories/psicometria"

# Print project directory to verify
print(f"Project directory: {project_directory}")
Project directory: /Users/corradocaudek/_repositories/psicometria

Capitolo 78

Esercizio 78.1

  1. Determinazione delle Matrici.

Le variabili indipendenti (\(x_1\), \(x_2\) e \(x_3\)) sono date, e possiamo creare la matrice delle variabili indipendenti \(\mathbf{X}\). Aggiungiamo una colonna di 1 per il termine costante \(\beta_0\):

\[ \mathbf{X} = \begin{pmatrix} 1 & 2 & 11 & 12 \\ 1 & 1 & 9 & 9 \\ 1 & 3 & 12 & 7 \\ 1 & 4 & 10 & 8 \\ 1 & 5 & 11 & 6 \end{pmatrix} \]

La matrice dei coefficienti \(\boldsymbol{\beta}\) è data come:

\[ \boldsymbol{\beta} = \begin{pmatrix} -1.402020 \\ 0.183838 \\ 1.405051 \\ -0.664646 \end{pmatrix} \]

La matrice degli errori \(\boldsymbol{\epsilon}\) non è fornita, ma la rappresentiamo comunque:

\[ \boldsymbol{\epsilon} = \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \end{pmatrix} \]

L’equazione del modello lineare in forma matriciale è:

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} \]

  1. Espansione del Modello.

Espandiamo il modello per ottenere le cinque equazioni esplicite utilizzando i valori forniti. Il modello per ogni osservazione \(i\) è:

\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \epsilon_i \]

Calcoliamo le cinque equazioni:

\[ \begin{aligned} y_1 &= -1.402020 + 0.183838 \cdot 2 + 1.405051 \cdot 11 - 0.664646 \cdot 12 + \epsilon_1 \\ y_2 &= -1.402020 + 0.183838 \cdot 1 + 1.405051 \cdot 9 - 0.664646 \cdot 9 + \epsilon_2 \\ y_3 &= -1.402020 + 0.183838 \cdot 3 + 1.405051 \cdot 12 - 0.664646 \cdot 7 + \epsilon_3 \\ y_4 &= -1.402020 + 0.183838 \cdot 4 + 1.405051 \cdot 10 - 0.664646 \cdot 8 + \epsilon_4 \\ y_5 &= -1.402020 + 0.183838 \cdot 5 + 1.405051 \cdot 11 - 0.664646 \cdot 6 + \epsilon_5 \\ \end{aligned} \]

  1. Calcolo dei Valori Predetti (\(\hat{y}_i\)).

Calcoliamo i valori predetti \(\hat{y}_i\) sostituendo \(\epsilon_i\) con 0:

\[ \begin{aligned} \hat{y}_1 &= -1.402020 + 0.183838 \cdot 2 + 1.405051 \cdot 11 - 0.664646 \cdot 12 \\ \hat{y}_2 &= -1.402020 + 0.183838 \cdot 1 + 1.405051 \cdot 9 - 0.664646 \cdot 9 \\ \hat{y}_3 &= -1.402020 + 0.183838 \cdot 3 + 1.405051 \cdot 12 - 0.664646 \cdot 7 \\ \hat{y}_4 &= -1.402020 + 0.183838 \cdot 4 + 1.405051 \cdot 10 - 0.664646 \cdot 8 \\ \hat{y}_5 &= -1.402020 + 0.183838 \cdot 5 + 1.405051 \cdot 11 - 0.664646 \cdot 6 \\ \end{aligned} \]

Calcoliamo i valori numerici:

\[ \begin{aligned} \hat{y}_1 &= -1.402020 + 0.367676 + 15.455561 - 7.975752 \approx 6.445465 \\ \hat{y}_2 &= -1.402020 + 0.183838 + 12.645459 - 5.981814 \approx 5.445463 \\ \hat{y}_3 &= -1.402020 + 0.551514 + 16.860612 - 4.652522 \approx 11.357584 \\ \hat{y}_4 &= -1.402020 + 0.735352 + 14.05051 - 5.317168 \approx 8.066674 \\ \hat{y}_5 &= -1.402020 + 0.91919 + 15.455561 - 3.987876 \approx 10.985855 \\ \end{aligned} \]

I valori predetti sono quindi:

\[ \begin{aligned} \hat{y}_1 &\approx 6.445 \\ \hat{y}_2 &\approx 5.445 \\ \hat{y}_3 &\approx 11.358 \\ \hat{y}_4 &\approx 8.067 \\ \hat{y}_5 &\approx 10.986 \\ \end{aligned} \]

  1. Calcolo degli errori casuali (\(\boldsymbol{\epsilon}\)).

Gli errori casuali (\(\epsilon_i\)) si ottengono sottraendo i valori predetti (\(\hat{y}_i\)) dai valori osservati (\(y_i\)). Per ciascun \(i\):

\[ \epsilon_i = y_i - \hat{y}_i \]

Calcoliamo gli errori casuali per ciascuna osservazione:

\[ \begin{aligned} \epsilon_1 &= 5.7 - 6.445 \approx -0.745 \\ \epsilon_2 &= 4.7 - 5.445 \approx -0.745 \\ \epsilon_3 &= 12.6 - 11.358 \approx 1.242 \\ \epsilon_4 &= 10.8 - 8.067 \approx 2.733 \\ \epsilon_5 &= 8.5 - 10.986 \approx -2.486 \\ \end{aligned} \]

Quindi, la matrice degli errori \(\boldsymbol{\epsilon}\) è:

\[ \boldsymbol{\epsilon} = \begin{pmatrix} -0.745 \\ -0.745 \\ 1.242 \\ 2.733 \\ -2.486 \end{pmatrix} \]

Esercizio 78.2

Ecco come implementare la formula dei minimi quadrati in Python utilizzando la libreria numpy e come verificare il risultato usando pingouin.linear_regression.

import numpy as np
import pandas as pd
import pingouin as pg

# Definizione della matrice X e del vettore y
X = np.array([[2, 11, 12], [1, 9, 9], [3, 12, 7], [4, 10, 8], [5, 11, 6]])
y = np.array([5.7, 4.7, 12.6, 10.8, 8.5])

# Implementazione con NumPy
X_numpy = np.hstack([np.ones((X.shape[0], 1)), X])
beta_numpy = np.linalg.inv(X_numpy.T @ X_numpy) @ X_numpy.T @ y
print("I coefficienti beta calcolati con NumPy sono:", beta_numpy)

# Implementazione con Pingouin
data = pd.DataFrame(X, columns=["x1", "x2", "x3"])
data["y"] = y
model_pingouin = pg.linear_regression(
    data[["x1", "x2", "x3"]], data["y"], add_intercept=True
)
print("\nI coefficienti beta calcolati con Pingouin sono:")
print(model_pingouin)

# Confronto dei risultati
print("\nConfrontiamo i coefficienti:")
print("NumPy:   ", beta_numpy)
print("Pingouin:", model_pingouin["coef"].values)
I coefficienti beta calcolati con NumPy sono: [-1.4020202   0.18383838  1.40505051 -0.66464646]

I coefficienti beta calcolati con Pingouin sono:
       names      coef         se         T      pval        r2   adj_r2  \
0  Intercept -1.402020  22.592392 -0.062057  0.960544  0.632638 -0.46945   
1         x1  0.183838   1.901446  0.096683  0.938640  0.632638 -0.46945   
2         x2  1.405051   1.960075  0.716835  0.604064  0.632638 -0.46945   
3         x3 -0.664646   1.214501 -0.547259  0.681221  0.632638 -0.46945   

     CI[2.5%]   CI[97.5%]  
0 -288.465573  285.661533  
1  -23.976322   24.343999  
2  -23.500063   26.310164  
3  -16.096345   14.767052  

Confrontiamo i coefficienti:
NumPy:    [-1.4020202   0.18383838  1.40505051 -0.66464646]
Pingouin: [-1.4020202   0.18383838  1.40505051 -0.66464646]

Capitolo 79

Esercizio 79.1

import os
import pandas as pd
import bambi as bmb
import arviz as az

# Define directories
home_directory = os.path.expanduser("~")
project_directory = f"{home_directory}/_repositories/psicometria"

data_path = os.path.join(project_directory, "data", "Aungle_Langer_2023.csv")
d = pd.read_csv(data_path)
d["Condition"] = d["Condition"].astype("category")
d.head()

# Define and fit the basic mixed-effects model
model = bmb.Model(
    "Healing ~ Condition + (1|Subject) + (1|ResponseId)",
    data=d
)

idata = model.fit(nuts_sampler="numpyro")
az.summary(idata, round_to=2, var_names="Condition")

model_varying_slopes = bmb.Model(
    "Healing ~ Condition + (1 + Condition|Subject) + (1 + Condition|ResponseId)",
    data=d,
)
idata_varying_slopes = model_varying_slopes.fit(nuts_sampler="numpyro")
az.summary(idata_varying_slopes, round_to=2, var_names="Condition")

Gelman & Brown (2024) hanno replicato l’analisi degli autori originali utilizzando un modello a effetti misti con pendenze casuali. Questa nuova analisi ha prodotto risultati sostanzialmente diversi. Nel nuovo modello statistico, i t-score sono stati notevolmente ridotti: da 10.7 a 3.0 per il confronto tra la condizione di 56 minuti e quella di 14 minuti, e da 2.5 a 0.7 per il confronto tra la condizione di 28 minuti e quella di 14 minuti. Inoltre, l’effetto della manipolazione sulla guarigione è stato stimato come altamente variabile tra i singoli soggetti, risultando a volte positivo e a volte negativo.

Gelman & Brown (2024) esprimono scetticismo riguardo alla capacità di questo studio di rivelare effetti reali del tempo percepito sulla guarigione fisica, basandosi su quattro argomentazioni principali.

  1. Selezione dei risultati: Il risultato statisticamente significativo emerso è solo uno dei molti confronti possibili. Lo studio ha raccolto dati su diverse variabili (ansia, stress, depressione, mindfulness, umore e tratti della personalità), aprendo la possibilità a numerose altre analisi. In assenza di una preregistrazione dello studio, non è possibile determinare quali analisi sarebbero state condotte se i dati fossero stati diversi.

  2. Variabilità tra i soggetti: La grande variazione stimata nell’effetto tra i soggetti suggerisce che qualsiasi effetto medio stimato sarebbe altamente dipendente dal campione specifico di partecipanti. Non vi è motivo di ritenere che i 33 partecipanti siano rappresentativi di una popolazione più ampia di interesse.

  3. Variabilità situazionale: Ci si aspetterebbe che l’effetto vari non solo tra le persone, ma anche tra le situazioni, sollevando preoccupazioni sulla stabilità e sulla generalizzabilità dei risultati. Inoltre, l’esperimento ha coinvolto solo contusioni molto lievi, non vere e proprie malattie o lesioni. Dunque, è difficile concludere, come fanno gli autori, che la tecnica del cupping abbia delle vere e propri proprietà curative.

  4. Mancanza di un meccanismo chiaro: Non è chiaro come la percezione del tempo possa influenzare la guarigione della pelle in questo contesto. Gli autori fanno riferimento a concetti generali come “l’unità mente-corpo”, ma non sono stati esaminati meccanismi specifici. Attribuire le differenze osservate al “tempo percepito” piuttosto che ad altri fattori che variavano tra le condizioni sembra essere una conclusione azzardata.

In conclusione, Gelman & Brown (2024) suggeriscono che, data la molteplicità di variabili in gioco nell’esperimento, potrebbe essere più plausibile interpretare i risultati come un esempio della capacità dei ricercatori di trovare pattern nel rumore in studi non adeguatamente controllati, piuttosto che come una prova dell’effetto del tempo percepito sulla guarigione fisica.