Obiettivo. Stimare l’effetto di class size (
stratio
) e % di English learners (english
) suscore
(media di matematica e lettura), controllando per omitted-variable bias. Mostriamo anche: calcolo manuale di OLS, \(R^2\) e \(R^2\) aggiustato, dummy trap, e le tre assunzioni chiave OLS con esempi simulati.
library(AER) # dataset CASchools
# library(car) # opzionale: per VIF
data("CASchools")
# Costruzione variabili chiave
CASchools$stratio <- with(CASchools, students/teachers) # proxy class size
CASchools$score <- with(CASchools, (math + read)/2) # outcome
head(CASchools[, c("district", "score", "stratio", "english")], 5)
Nota didattica. stratio
misura il
rapporto studenti/docenti (più alto = classi più grandi).
score
è la media dei punteggi di matematica e lettura.
Commento: Le variabili principali descrivono
l’ambiente scolastico: stratio
indica la dimensione media
della classe (più alta = classi più grandi), mentre score
misura la performance scolastica media. La percentuale di studenti non
madrelingua (english
) sarà un potenziale fattore
esplicativo della performance.
# (Opzionale) scatter matrix per ispezione pattern e outlier
pairs(CASchools[, c("score", "stratio", "english")])
# Correlazioni grezze (indicative, non causali)
cor(CASchools$stratio, CASchools$english)
## [1] 0.1876424
cor(CASchools$score, CASchools$english)
## [1] -0.6441238
Lettura. Di norma si osserva: (i) score
e english
negativamente correlati; (ii) score
e stratio
negativamente correlati; (iii)
stratio
e english
positivamente correlati
(potenziale confondente).
Commento: La correlazione positiva tra
stratio
e english
(≈ 0,19) suggerisce che le
scuole con classi più grandi hanno più studenti non madrelingua.
score
e english
mostrano invece una
correlazione negativa (≈ −0,64): dove ci sono più studenti con
difficoltà linguistiche, i punteggi medi sono inferiori.
reg_mult <- lm(score ~ stratio + english, data = CASchools)
summary(reg_mult)
##
## Call:
## lm(formula = score ~ stratio + english, data = CASchools)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.845 -10.240 -0.308 9.815 43.461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 686.03224 7.41131 92.566 < 2e-16 ***
## stratio -1.10130 0.38028 -2.896 0.00398 **
## english -0.64978 0.03934 -16.516 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared: 0.4264, Adjusted R-squared: 0.4237
## F-statistic: 155 on 2 and 417 DF, p-value: < 2.2e-16
Interpretazione. Il coefficiente su
stratio
misura la variazione attesa di score
al crescere della dimensione della classe a parità di
english
. Il coefficiente su english
misura
l’effetto della quota di English learners a parità di
dimensione della classe.
Suggerimento: per diagnosi di multicollinearità, è possibile usare
car::vif(reg_mult)
(facoltativo).
reg_semp <- lm(score ~ stratio, data = CASchools)
summary(reg_semp)
##
## Call:
## lm(formula = score ~ stratio, data = CASchools)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.727 -14.251 0.483 12.822 48.540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 698.9329 9.4675 73.825 < 2e-16 ***
## stratio -2.2798 0.4798 -4.751 2.78e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared: 0.05124, Adjusted R-squared: 0.04897
## F-statistic: 22.58 on 1 and 418 DF, p-value: 2.783e-06
Commento: * Nel modello semplice
(score ~ stratio
), l’effetto di stratio
è
−2,28 (p < 0,001), ma il modello spiega solo il 5% della variabilità
del punteggio (R² = 0,05). * Nel modello multiplo
(score ~ stratio + english
), i coefficienti diventano:
stratio
−1,10 (p = 0,004) e english
−0,65 (p
< 0,001). La varianza spiegata sale al 42,6%. * Questo mostra un
omitted variable bias nel modello semplice: parte
dell’effetto apparente di stratio
era in realtà dovuto alla
quota di studenti EL. * L’aggiunta di english
migliora la
capacità esplicativa del modello senza introdurre collinearità (VIF ≈
1,03).
# OLS: beta = (X'X)^{-1} X'Y
X <- model.matrix(~ stratio + english, data = CASchools) # include intercetta
Y <- CASchools$score
XtX <- t(X) %*% X
XtY <- t(X) %*% Y
beta_hat <- solve(XtX, XtY)
beta_hat
## [,1]
## (Intercept) 686.0322445
## stratio -1.1012956
## english -0.6497768
coef(reg_mult) # confronto con lm()
## (Intercept) stratio english
## 686.0322445 -1.1012956 -0.6497768
Spiegazione. model.matrix()
genera la
matrice del modello con intercetta. Confrontiamo il risultato con
coef(lm(...))
per verificare l’identità con la formula
OLS.
Commento: I risultati coincidono con
lm()
, confermando la formula OLS ( \(\hat{\beta} = (X'X)^{-1} X'Y\) ).
L’approccio matriciale è utile per visualizzare il meccanismo interno
della stima.
# già stimato come reg_semp sopra; qui solo richiamo al confronto numerico
summary(reg_semp)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 698.932949 9.4674911 73.824516 6.569846e-242
## stratio -2.279808 0.4798255 -4.751327 2.783308e-06
summary(reg_mult)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 686.0322445 7.41131160 92.565565 3.871327e-280
## stratio -1.1012956 0.38027827 -2.896026 3.978059e-03
## english -0.6497768 0.03934254 -16.515882 1.657448e-47
Nota. Il confronto tra semplice e multipla evidenzia
l’effetto del controllo per english
(riduzione del bias da
variabile omessa).
# Modello semplice
TSS_sempl <- sum( (CASchools$score - mean(CASchools$score))^2 )
SSR_sempl <- sum( residuals(reg_semp)^2 )
R2_sempl <- 1 - SSR_sempl/TSS_sempl
n <- nrow(CASchools); k_sempl <- length(coef(reg_semp))-1
R2adj_sempl <- 1 - ((n-1)/(n-k_sempl-1)) * SSR_sempl/TSS_sempl
c(R2_sempl = R2_sempl, R2adj_sempl = R2adj_sempl,
R2_sempl_summary = summary(reg_semp)$r.squared,
R2adj_sempl_summary = summary(reg_semp)$adj.r.squared)
## R2_sempl R2adj_sempl R2_sempl_summary R2adj_sempl_summary
## 0.05124009 0.04897033 0.05124009 0.04897033
# Modello multiplo
TSS_mult <- sum( (CASchools$score - mean(CASchools$score))^2 )
SSR_mult <- sum( residuals(reg_mult)^2 )
R2_mult <- 1 - SSR_mult/TSS_mult
k_mult <- length(coef(reg_mult))-1
R2adj_mult <- 1 - ((n-1)/(n-k_mult-1)) * SSR_mult/TSS_mult
c(R2_mult = R2_mult, R2adj_mult = R2adj_mult,
R2_mult_summary = summary(reg_mult)$r.squared,
R2adj_mult_summary = summary(reg_mult)$adj.r.squared)
## R2_mult R2adj_mult R2_mult_summary R2adj_mult_summary
## 0.4264315 0.4236805 0.4264315 0.4236805
Lettura. \(R^2\)
cresce aggiungendo regressori; \(R^2\)
aggiustato penalizza regressori non informativi. Usiamo
summary()
come controllo.
par(mfrow=c(1,2))
hist(CASchools$income, main="income (raw)")
hist(log(CASchools$income), main="log(income)")
par(mfrow=c(1,1))
Nota. income
è spesso asimmetrica;
log(income)
può migliorare linearità e
omoschedasticità.
# Modello corretto con county come fattore (usa K-1 dummies)
reg_mult_dummy_ok <- lm(score ~ stratio + english + county, data = CASchools)
summary(reg_mult_dummy_ok)[1:5]
## $call
## lm(formula = score ~ stratio + english + county, data = CASchools)
##
## $terms
## score ~ stratio + english + county
## attr(,"variables")
## list(score, stratio, english, county)
## attr(,"factors")
## stratio english county
## score 0 0 0
## stratio 1 0 0
## english 0 1 0
## county 0 0 1
## attr(,"term.labels")
## [1] "stratio" "english" "county"
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(score, stratio, english, county)
## attr(,"dataClasses")
## score stratio english county
## "numeric" "numeric" "numeric" "factor"
##
## $residuals
## 1 2 3 4 5
## 2.184919e-13 1.484800e+01 1.096092e+01 -6.585082e+00 -2.624785e+00
## 6 7 8 9 10
## -3.529752e+01 -3.206204e+00 -1.295080e+01 -1.822225e+01 -9.594270e+00
## 11 12 13 14 15
## -9.740427e+00 3.227147e+00 -8.686438e+00 -2.422734e+01 1.349916e+01
## 16 17 18 19 20
## -3.561364e+00 1.735263e+01 -7.711425e-01 1.064360e+01 -8.875678e-01
## 21 22 23 24 25
## -1.689278e+01 -1.897349e+01 -3.150168e+00 -1.769754e-01 1.235101e+01
## 26 27 28 29 30
## -5.615932e+00 1.240741e+01 1.503014e+01 5.902810e+00 -2.209054e+01
## 31 32 33 34 35
## 4.028713e+00 4.999102e+00 6.824888e+00 -1.967623e+00 2.190453e+00
## 36 37 38 39 40
## -1.114207e+01 1.635394e-01 -6.413698e+00 -7.039053e+00 8.709780e+00
## 41 42 43 44 45
## -2.978897e+00 -3.242228e+00 -5.791635e+00 1.293393e+01 -4.880403e+00
## 46 47 48 49 50
## -9.933843e+00 -1.400279e+01 -1.380067e+01 -5.704611e+00 -1.098121e+01
## 51 52 53 54 55
## -3.181539e+00 -3.105056e+00 -6.436673e+00 -1.592219e+01 -8.254160e+00
## 56 57 58 59 60
## -1.625606e+01 -7.850513e+00 4.388328e+00 -5.961293e+00 -7.725731e+00
## 61 62 63 64 65
## 2.777986e+00 2.267905e+00 5.433906e-02 -9.987874e+00 -1.238481e+01
## 66 67 68 69 70
## -1.773346e+01 -8.129408e+00 -8.815566e+00 2.173657e-03 -7.653437e-01
## 71 72 73 74 75
## -1.475280e+01 -1.845516e+01 -1.118936e+01 -3.592364e-01 -4.339916e+00
## 76 77 78 79 80
## -8.005128e+00 -4.033878e+01 -2.063854e+01 -1.705657e+01 -1.839235e+01
## 81 82 83 84 85
## -8.294629e+00 -1.061065e+01 -9.747193e-01 -5.992079e+00 2.443730e+00
## 86 87 88 89 90
## -2.166191e+01 8.381554e-01 1.405878e+01 -8.173860e+00 -8.149898e+00
## 91 92 93 94 95
## -1.458526e+01 -3.337281e+00 -1.717233e+01 -1.529281e+01 -2.891315e+01
## 96 97 98 99 100
## -9.794971e+00 -1.340420e+01 5.864899e-01 4.669948e+00 -5.490616e+00
## 101 102 103 104 105
## -7.806241e+00 -7.859820e+00 -1.235999e+01 4.336626e-15 -7.559900e+00
## 106 107 108 109 110
## -5.239048e+00 -8.355907e+00 -1.306238e+00 -1.421131e+01 -2.008488e+01
## 111 112 113 114 115
## 1.507931e+01 -1.582352e+01 -9.970709e+00 -4.909598e+00 -7.978804e+00
## 116 117 118 119 120
## -1.087845e+01 1.614604e+01 -5.979303e+00 -1.626854e+01 8.689882e+00
## 121 122 123 124 125
## -1.181248e+01 -2.715055e+00 3.332791e+00 -1.912961e+01 -5.459091e+00
## 126 127 128 129 130
## -1.449629e-01 3.027267e-01 1.129168e+01 -2.012869e+00 -2.156641e+01
## 131 132 133 134 135
## -8.204846e+00 -2.846798e+00 -1.835015e+00 -1.069686e+01 -2.042341e+00
## 136 137 138 139 140
## -1.684108e+01 -7.961356e+00 -6.410096e+00 6.259687e+00 1.671526e+00
## 141 142 143 144 145
## -9.632995e-01 -7.293584e+00 8.193245e+00 -1.492108e+01 -1.590432e+01
## 146 147 148 149 150
## -1.860278e+01 3.337281e+00 -1.432012e+01 1.043318e+01 -8.584065e+00
## 151 152 153 154 155
## -2.700001e+00 -1.096992e+00 -6.163047e+00 -6.128183e+00 -6.826072e-01
## 156 157 158 159 160
## -2.199789e+00 2.869712e+00 2.628800e+00 2.943810e+00 -4.486631e+00
## 161 162 163 164 165
## -1.406387e+01 -8.086376e+00 -2.772672e+00 -3.256685e-01 4.607539e-01
## 166 167 168 169 170
## -2.858299e+00 5.675126e-02 -3.194580e+00 -1.819518e-01 -4.232669e+00
## 171 172 173 174 175
## -1.412591e+01 -4.438665e+00 -3.925993e+00 -6.581436e+00 3.510317e+00
## 176 177 178 179 180
## -5.489307e+00 -1.523441e+01 -5.552881e+00 -1.698283e+01 -2.394406e+01
## 181 182 183 184 185
## -2.039799e+00 -2.321471e-01 -9.316364e+00 -4.532274e+00 3.363327e-02
## 186 187 188 189 190
## -8.840568e+00 -2.104286e+01 -1.128911e+01 -5.712246e+00 -5.840865e+00
## 191 192 193 194 195
## -2.846221e-01 -1.097358e+01 3.588466e-01 -5.794477e+00 3.141566e+00
## 196 197 198 199 200
## -7.064020e+00 -5.006274e+00 -1.754128e+01 -1.491602e+01 4.187211e+00
## 201 202 203 204 205
## 5.027730e-01 -7.815772e+00 1.202550e+00 2.606968e+00 2.542033e+00
## 206 207 208 209 210
## -6.043178e-01 -8.003840e+00 1.324479e+01 -1.069581e+01 1.310776e+01
## 211 212 213 214 215
## -2.215211e+00 1.855177e+01 -1.510216e+00 -1.449364e+01 2.697611e+00
## 216 217 218 219 220
## 1.298134e+01 -1.697250e+01 5.692109e+00 1.938535e+00 -2.917344e-01
## 221 222 223 224 225
## 5.165410e+00 -6.354947e+00 6.167789e+00 -8.312541e+00 6.749264e+00
## 226 227 228 229 230
## -6.121746e+00 2.940023e+00 -2.102116e+00 -1.500087e+01 -1.776981e+01
## 231 232 233 234 235
## 4.353513e+00 -2.310923e-01 5.391338e-15 -3.351222e-01 -5.894547e+00
## 236 237 238 239 240
## -7.625666e+00 9.236735e+00 -7.343697e+00 -1.665511e+00 -1.192123e+01
## 241 242 243 244 245
## -1.180259e+01 2.027767e+00 5.256390e-01 4.510684e+00 1.989037e+00
## 246 247 248 249 250
## -1.080560e-01 6.840197e+00 -6.345772e+00 -6.040299e+00 -5.399156e+00
## 251 252 253 254 255
## 7.516380e+00 3.059870e-15 -8.214217e+00 -3.385728e-03 1.155493e+01
## 256 257 258 259 260
## 4.998395e+00 -8.607192e+00 1.485904e+01 8.705424e+00 4.383354e+00
## 261 262 263 264 265
## -2.971564e-01 3.113546e+01 -3.391654e-01 3.860093e+00 -5.452313e+00
## 266 267 268 269 270
## 1.059197e+01 -1.043734e+01 7.869980e+00 1.201799e+00 1.803732e+00
## 271 272 273 274 275
## -2.901994e+00 4.712967e+00 1.439349e+01 -8.535954e-01 6.326991e+00
## 276 277 278 279 280
## 6.226195e+00 -6.696035e+00 9.390111e+00 3.515851e+00 -5.561145e+00
## 281 282 283 284 285
## 1.370427e+01 -3.349940e+00 2.800912e+00 4.847172e+00 4.215492e-01
## 286 287 288 289 290
## 6.176116e+00 5.285714e+00 1.884603e+00 1.148650e+00 -1.797783e+00
## 291 292 293 294 295
## -1.648920e+00 7.594294e+00 1.207642e+01 -1.012306e+00 -2.113288e-01
## 296 297 298 299 300
## 2.986747e+00 7.194675e+00 -1.304736e+00 3.538896e+00 8.851116e+00
## 301 302 303 304 305
## 3.325633e+00 6.214208e+00 -7.698582e+00 5.237191e-01 -1.796977e+00
## 306 307 308 309 310
## -5.890219e+00 -6.497326e+00 1.645893e+00 -3.054585e+00 2.391130e+00
## 311 312 313 314 315
## 1.419691e+01 -1.178672e+00 9.553984e+00 3.990792e+00 1.057784e+01
## 316 317 318 319 320
## 1.535518e+01 -2.373313e+00 -5.986953e+00 -3.751816e-01 9.604926e+00
## 321 322 323 324 325
## 1.397130e+01 1.119546e+01 -4.605756e-01 2.236685e+01 1.603298e+01
## 326 327 328 329 330
## -6.078500e+00 4.776655e+00 2.132144e+01 1.212568e+01 2.275511e+00
## 331 332 333 334 335
## 1.157663e+01 1.625793e+01 -1.510344e+01 -8.256890e+00 1.647989e+01
## 336 337 338 339 340
## 1.198238e+01 -9.698789e-01 1.208335e+01 1.567119e+01 -6.716096e+00
## 341 342 343 344 345
## 4.711708e+00 5.013637e+00 6.728280e+00 6.522520e+00 5.662579e+00
## 346 347 348 349 350
## 1.357917e+01 1.144410e+01 5.958386e+00 5.423207e+00 3.568553e+00
## 351 352 353 354 355
## -3.553559e+00 5.975067e+00 1.027824e+01 6.345772e+00 1.047799e+01
## 356 357 358 359 360
## 2.790272e+01 6.856566e+00 7.610472e+00 -3.481053e+00 7.449980e+00
## 361 362 363 364 365
## 5.104245e+00 9.482730e+00 3.297287e+00 8.103170e+00 4.605209e+00
## 366 367 368 369 370
## 1.709121e+01 3.609716e-01 1.229331e+01 9.456444e+00 2.096897e+01
## 371 372 373 374 375
## -3.609716e-01 1.784951e+01 1.288489e+01 1.927885e+01 2.918300e+00
## 376 377 378 379 380
## 1.388293e+01 -4.415097e+00 -1.386952e+00 1.589297e+01 1.112850e+01
## 381 382 383 384 385
## 1.608869e+01 1.170963e+01 4.636373e-02 1.324867e+01 1.274613e+01
## 386 387 388 389 390
## 1.370986e+01 -2.462048e+00 4.381598e+00 1.388600e+01 -2.659429e-01
## 391 392 393 394 395
## 3.028612e+00 1.663386e+01 2.332683e+01 4.852342e+00 3.139115e+00
## 396 397 398 399 400
## 1.341766e+01 1.207315e+01 3.167283e+00 1.510207e+01 1.204219e+01
## 401 402 403 404 405
## 1.054999e+01 6.859179e+00 2.017582e+01 1.572856e+01 8.646676e+00
## 406 407 408 409 410
## 1.076769e+01 1.927958e+01 1.973401e+01 1.426620e+01 2.060571e+01
## 411 412 413 414 415
## 1.924274e+01 1.727142e+01 1.485346e+01 1.733468e+01 2.122927e+01
## 416 417 418 419 420
## 2.701450e+01 2.829978e+01 -4.502052e+00 8.222972e+00 -8.222972e+00
##
## $coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 711.6842038 13.03505503 54.5977138 5.618271e-180
## stratio -1.1673741 0.36871580 -3.1660538 1.672347e-03
## english -0.6693875 0.03630142 -18.4397069 1.973521e-54
## countyButte -37.1368313 12.15671248 -3.0548416 2.413671e-03
## countyCalaveras -29.4941873 15.94346921 -1.8499228 6.511542e-02
## countyContra Costa -12.0677199 12.03366140 -1.0028303 3.165929e-01
## countyEl Dorado -22.8805038 11.80464022 -1.9382635 5.334530e-02
## countyFresno -37.5413002 11.73937207 -3.1978968 1.502637e-03
## countyGlenn -19.6594048 13.05499761 -1.5058911 1.329415e-01
## countyHumboldt -31.0497435 11.56977732 -2.6836941 7.605787e-03
## countyImperial -25.6576608 12.25916276 -2.0929374 3.702986e-02
## countyInyo -22.3843871 15.93451602 -1.4047736 1.609208e-01
## countyKern -33.9080709 11.50458091 -2.9473539 3.406600e-03
## countyKings -31.2349076 11.88855116 -2.6273099 8.961344e-03
## countyLake -47.2031767 13.77225249 -3.4274115 6.773068e-04
## countyLassen -24.0299488 12.37180890 -1.9423149 5.285137e-02
## countyLos Angeles -24.4769938 11.55527183 -2.1182534 3.481440e-02
## countyMadera -28.5121217 12.33429026 -2.3116143 2.134321e-02
## countyMarin -5.8796729 11.92742813 -0.4929540 6.223352e-01
## countyMendocino -38.2846570 15.93606386 -2.4023910 1.677632e-02
## countyMerced -25.9804684 11.85338094 -2.1918192 2.901017e-02
## countyMonterey -24.8749185 12.14468620 -2.0482142 4.123864e-02
## countyNevada -23.6045613 11.85356393 -1.9913472 4.717170e-02
## countyOrange -13.9980309 11.85025934 -1.1812426 2.382590e-01
## countyPlacer -21.1644781 11.77525618 -1.7973688 7.308582e-02
## countyRiverside -30.9276995 12.65173894 -2.4445414 1.496581e-02
## countySacramento -38.1918441 12.05183573 -3.1689649 1.656127e-03
## countySan Benito -29.1147507 12.99765185 -2.2400008 2.567928e-02
## countySan Bernardino -29.6823023 11.85642483 -2.5034783 1.272449e-02
## countySan Diego -13.2746563 11.55783693 -1.1485416 2.514811e-01
## countySan Joaquin -32.9652669 12.17587106 -2.7074258 7.092668e-03
## countySan Luis Obispo -19.6674749 13.78043703 -1.4272026 1.543581e-01
## countySan Mateo -11.1536007 11.58303412 -0.9629257 3.362087e-01
## countySanta Barbara -6.3458247 11.79507437 -0.5380063 5.908935e-01
## countySanta Clara -9.2180287 11.56849987 -0.7968214 4.260616e-01
## countySanta Cruz -8.4160142 12.03515811 -0.6992857 4.848093e-01
## countyShasta -29.6967367 11.67641323 -2.5433098 1.138379e-02
## countySiskiyou -41.3486401 11.87658650 -3.4815256 5.576430e-04
## countySonoma -17.6618621 11.44515366 -1.5431739 1.236365e-01
## countyStanislaus -23.1836391 12.05854241 -1.9225905 5.529273e-02
## countySutter -34.7840243 12.18286502 -2.8551596 4.541959e-03
## countyTehama -30.7200476 11.94333854 -2.5721491 1.049319e-02
## countyTrinity -13.8183108 13.77029322 -1.0034870 3.162765e-01
## countyTulare -34.0546777 11.54639264 -2.9493781 3.384876e-03
## countyTuolumne -34.4739143 12.14373351 -2.8388234 4.775822e-03
## countyVentura -20.3917949 11.97740570 -1.7025219 8.949071e-02
## countyYuba -22.1379269 13.78518172 -1.6059220 1.091373e-01
##
## $aliased
## (Intercept) stratio english
## FALSE FALSE FALSE
## countyButte countyCalaveras countyContra Costa
## FALSE FALSE FALSE
## countyEl Dorado countyFresno countyGlenn
## FALSE FALSE FALSE
## countyHumboldt countyImperial countyInyo
## FALSE FALSE FALSE
## countyKern countyKings countyLake
## FALSE FALSE FALSE
## countyLassen countyLos Angeles countyMadera
## FALSE FALSE FALSE
## countyMarin countyMendocino countyMerced
## FALSE FALSE FALSE
## countyMonterey countyNevada countyOrange
## FALSE FALSE FALSE
## countyPlacer countyRiverside countySacramento
## FALSE FALSE FALSE
## countySan Benito countySan Bernardino countySan Diego
## FALSE FALSE FALSE
## countySan Joaquin countySan Luis Obispo countySan Mateo
## FALSE FALSE FALSE
## countySanta Barbara countySanta Clara countySanta Cruz
## FALSE FALSE FALSE
## countyShasta countySiskiyou countySonoma
## FALSE FALSE FALSE
## countyStanislaus countySutter countyTehama
## FALSE FALSE FALSE
## countyTrinity countyTulare countyTuolumne
## FALSE FALSE FALSE
## countyVentura countyYuba
## FALSE FALSE
# Creiamo una dummy "di troppo" (Alameda) -> collinearità perfetta
CASchools$Alameda_dummy <- ifelse(CASchools$county == "Alameda", 1, 0)
reg_mult_dummy_bad <- lm(score ~ stratio + english + county + Alameda_dummy,
data = CASchools)
summary(reg_mult_dummy_bad) # attesi NA / coeff. non stimabili
##
## Call:
## lm(formula = score ~ stratio + english + county + Alameda_dummy,
## data = CASchools)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.339 -6.613 -0.249 6.768 31.135
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 711.6842 13.0351 54.598 < 2e-16 ***
## stratio -1.1674 0.3687 -3.166 0.001672 **
## english -0.6694 0.0363 -18.440 < 2e-16 ***
## countyButte -37.1368 12.1567 -3.055 0.002414 **
## countyCalaveras -29.4942 15.9435 -1.850 0.065115 .
## countyContra Costa -12.0677 12.0337 -1.003 0.316593
## countyEl Dorado -22.8805 11.8046 -1.938 0.053345 .
## countyFresno -37.5413 11.7394 -3.198 0.001503 **
## countyGlenn -19.6594 13.0550 -1.506 0.132941
## countyHumboldt -31.0497 11.5698 -2.684 0.007606 **
## countyImperial -25.6577 12.2592 -2.093 0.037030 *
## countyInyo -22.3844 15.9345 -1.405 0.160921
## countyKern -33.9081 11.5046 -2.947 0.003407 **
## countyKings -31.2349 11.8886 -2.627 0.008961 **
## countyLake -47.2032 13.7722 -3.427 0.000677 ***
## countyLassen -24.0299 12.3718 -1.942 0.052851 .
## countyLos Angeles -24.4770 11.5553 -2.118 0.034814 *
## countyMadera -28.5121 12.3343 -2.312 0.021343 *
## countyMarin -5.8797 11.9274 -0.493 0.622335
## countyMendocino -38.2847 15.9361 -2.402 0.016776 *
## countyMerced -25.9805 11.8534 -2.192 0.029010 *
## countyMonterey -24.8749 12.1447 -2.048 0.041239 *
## countyNevada -23.6046 11.8536 -1.991 0.047172 *
## countyOrange -13.9980 11.8503 -1.181 0.238259
## countyPlacer -21.1645 11.7753 -1.797 0.073086 .
## countyRiverside -30.9277 12.6517 -2.445 0.014966 *
## countySacramento -38.1918 12.0518 -3.169 0.001656 **
## countySan Benito -29.1148 12.9977 -2.240 0.025679 *
## countySan Bernardino -29.6823 11.8564 -2.503 0.012724 *
## countySan Diego -13.2747 11.5578 -1.149 0.251481
## countySan Joaquin -32.9653 12.1759 -2.707 0.007093 **
## countySan Luis Obispo -19.6675 13.7804 -1.427 0.154358
## countySan Mateo -11.1536 11.5830 -0.963 0.336209
## countySanta Barbara -6.3458 11.7951 -0.538 0.590894
## countySanta Clara -9.2180 11.5685 -0.797 0.426062
## countySanta Cruz -8.4160 12.0352 -0.699 0.484809
## countyShasta -29.6967 11.6764 -2.543 0.011384 *
## countySiskiyou -41.3486 11.8766 -3.482 0.000558 ***
## countySonoma -17.6619 11.4451 -1.543 0.123637
## countyStanislaus -23.1836 12.0585 -1.923 0.055293 .
## countySutter -34.7840 12.1829 -2.855 0.004542 **
## countyTehama -30.7201 11.9433 -2.572 0.010493 *
## countyTrinity -13.8183 13.7703 -1.003 0.316276
## countyTulare -34.0547 11.5464 -2.949 0.003385 **
## countyTuolumne -34.4739 12.1437 -2.839 0.004776 **
## countyVentura -20.3918 11.9774 -1.703 0.089491 .
## countyYuba -22.1379 13.7852 -1.606 0.109137
## Alameda_dummy NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.24 on 373 degrees of freedom
## Multiple R-squared: 0.69, Adjusted R-squared: 0.6518
## F-statistic: 18.05 on 46 and 373 DF, p-value: < 2.2e-16
Perché succede? Con \(K\) categorie servono \(K-1\) dummies. Aggiungere la K-esima rende una colonna combinazione lineare delle altre (collinearità perfetta).
Commento: Inserendo le dummies per
county
, il modello spiega il 69% della varianza totale (R²
= 0,69). stratio
(−1,17) e english
(−0,67)
restano negativi e significativi anche controllando per eterogeneità
territoriale. L’aggiunta della dummy Alameda_dummy
genera
collinearità perfetta: R la segnala con un coefficiente NA, dimostrando
la dummy trap (una dummy per ogni categoria + intercetta =
colonna ridondante).
set.seed(123)
N <- 100
X <- runif(N, -5, 5)
u <- rnorm(N, 1)
Y <- 1 + X^2 + 2*X + u # verità con termine quadratico
reg_lin <- lm(Y ~ X) # modello mal specificato
summary(reg_lin)
##
## Call:
## lm(formula = Y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.475 -6.461 -1.244 6.027 16.988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.9883 0.7321 13.642 < 2e-16 ***
## X 2.0590 0.2582 7.975 2.87e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.321 on 98 degrees of freedom
## Multiple R-squared: 0.3936, Adjusted R-squared: 0.3874
## F-statistic: 63.6 on 1 and 98 DF, p-value: 2.871e-12
prediction <- predict(lm(Y ~ X + I(X^2)), data.frame(X = sort(X)))
plot(Y ~ X, main="Violazione E(u|X)=0: modello lineare vs vero (quadratico)")
abline(reg_lin, col="red", lwd=2)
lines(sort(X), prediction, col="blue", lwd=2)
legend("topleft", c("Lineare (stimato)","Vero (X + X^2)"),
col=c("red","blue"), lty=1, bty="n")
plot(X, resid(reg_lin), main="Residui vs X (struttura non casuale)")
abline(h=0)
lines(lowess(X, resid(reg_lin)), col="red", lwd=2)
mean(resid(reg_lin)) # in media zero, ma non condizionatamente a X
## [1] 3.017551e-17
Conclusione. Modello mal specificato → residui
strutturati → coefficiente su X
distorto.
N <- 1000; set.seed(1)
X <- runif(N, -5, 5); u <- rnorm(N, 1)
Y <- 1 + X^2 + 2*X + u
db <- data.frame(Y, X, X2 = X^2)
R <- 2000; N_sub <- 25 # (ridotto a 2000 per tempi di knitting)
beta_1_hat <- rep(NA, R)
for(i in 1:R){
idx <- sample(1:N, N_sub, replace=TRUE)
reg <- lm(Y ~ X, data=db[idx,])
beta_1_hat[i] <- coef(reg)[2]
}
mean(beta_1_hat - 2) # bias rispetto al vero coefficiente su X (2)
## [1] 0.1127015
N <- 1000; set.seed(2)
X1 <- runif(N, -5, 5)
X2 <- X1 + rnorm(N, sd=10) # correlazione imperfetta tra regressori
cor(X1, X2)
## [1] 0.3103858
u <- rnorm(N, 1)
Y <- 1 + 0.5*X1 - 4*X2 + u
db <- data.frame(Y, X1, X2)
R <- 2000; N_sub <- 25 # (ridotto a 2000 per tempi)
beta_1_hat <- rep(NA, R)
for(i in 1:R){
idx <- sample(1:N, N_sub, replace=TRUE)
reg <- lm(Y ~ X1, data=db[idx,])
beta_1_hat[i] <- coef(reg)[2]
}
mean(beta_1_hat - 0.5) # bias su beta1 quando si omette X2
## [1] -4.393373
Date <- seq(as.Date("1971/1/1"), as.Date("2020/1/1"), by="years")
X <- c(5000, rep(NA, length(Date)-1))
set.seed(10)
for(i in 2:length(Date)){
X[i] <- -50 + 0.98*X[i-1] + rnorm(1, sd=200) # dipendenza seriale
}
plot(x=Date, y=X, type="l", col="steelblue",
ylab="X", xlab="Time",
main="Esempio di autocorrelazione (non-i.i.d.)")
X <- sort(runif(50, 30, 70))
set.seed(3)
Y <- rnorm(50, mean=200, sd=50)
Y[40] <- 2000 # outlier artificiale
fit <- lm(Y ~ X)
fit_no_outlier <- lm(Y[-40] ~ X[-40])
plot(Y ~ X, ylim=c(0,2000), main="Effetto di un outlier sulla stima OLS")
abline(fit, lwd=2)
abline(fit_no_outlier, col="red", lwd=2)
legend("topleft", c("Con outlier","Senza outlier"),
col=c("black","red"), lty=1, bty="n")
Commento: * Assunzione 1 – E(u|X)=0: stimando un modello lineare su una relazione quadratica, i residui non sono casuali rispetto a X → violazione di specificazione. La simulazione mostra un bias medio su β̂_X ≈ +0,10. * Assunzione 1 – Omitted variable bias: quando Y dipende da X1 e X2, ma si stima solo su X1, il bias medio su β̂_X1 è −4,35: la stima è gravemente distorta. * Assunzione 2 – i.i.d.: la simulazione in serie storiche mostra autocorrelazione; la dipendenza seriale viola l’indipendenza richiesta da OLS. * Assunzione 3 – Outlier: un singolo outlier altera significativamente la retta stimata, evidenziando che OLS non è robusto.
stratio
controllando per english
,
mitigando il bias da variabile omessa.Risultati chiave: * stratio
e
english
influenzano negativamente score
;
controllare per english
dimezza l’effetto di
stratio
(da −2,28 a −1,10). * Le fisse di contea
migliorano notevolmente la capacità esplicativa del modello (R² da 0,42
a 0,69). * Nessuna multicollinearità significativa (VIF ≈ 1).
Implicazioni metodologiche: * L’inclusione di
variabili di controllo pertinenti (come english
) è cruciale
per evitare bias. * Verificare sempre le assunzioni OLS: specificazione
corretta, indipendenza e assenza di outlier. * Le simulazioni dimostrano
l’importanza della corretta specificazione funzionale e della robustezza
delle stime.
Conclusione: Il dataset CASchools offre un chiaro esempio di come la regressione multipla permetta di distinguere effetti reali da correlazioni spurie. Una buona pratica econometrica richiede sempre il controllo per variabili rilevanti e la verifica delle ipotesi del modello per garantire inferenze affidabili.