Zobecněné lineární modely
Problém 1: Ceny nemovitostí Poznámky k řešení1
Zadání 1. Mají nemovitosti určené k bydlení vyšší cenu tam, kde je čistší ovzduší? Pokud ano, o kolik? 2. Lze vztah mezi znečištěním a cenou, pokud existuje, vysvětlit tím, že ve znečištěných oblastech bydlí chudší lidé, menšiny, jsou tam horší veřejné služby, atd.? 3. Myslíte, že cílený program na zlepšení čistoty ovzduší by vedl ke zvýšení cen rodinných domků v dané lokalitě?
Postup • Načtu data: data1 <- read.csv("cvic1.csv") Ověřím si velikost dat a jména veličin: names(data1); dim(data1) Vypíšu si základní popisné charakteristiky veličin: summary(data1) Vidím, že (i) v datech nejsou chybějící hodnoty; (ii) všechny veličiny jsou spojité kromě chas, která je nula-jedničková. Zajistím si přímý přístup k veličinám: attach(data1). • Podívám se na nejdůležitější veličiny podrobněji. Např. histogramy (hist(medv), hist(nox)), tabulky četností: > table(cut(nox,c(-Inf,seq(0.4,0.8,by=0.1),Inf))) (-Inf,0.4] (0.4,0.5] (0.5,0.6] (0.6,0.7] (0.7,0.8] 11 181 149 104 45 > table(cut(medv,c(-Inf,seq(5,30,by=5),Inf))) (-Inf,5] 2
(5,10] 22
(10,15] 73
(15,20] 118
(20,25] 167
(0.8,Inf] 16
(25,30] (30,Inf] 40 84
Prozkoumám popisně vztah mezi nox a medv. Např. obrázek (scatterplot ) vyhlazený neparametrickou křivkou lowess: (plot(nox,medv); lines(lowess(nox,medv))) nebo tabulku průměrů medv podle intervalů nox: > tapply(medv,cut(nox,c(-Inf,seq(0.4,0.8,by=0.1),Inf)),mean) (-Inf,0.4] (0.4,0.5] (0.5,0.6] (0.6,0.7] (0.7,0.8] (0.8,Inf] 25.48182 26.68729 21.70268 19.92885 16.04000 16.42500 • Uvědomím si, že nox nabývá hodnot zhruba mezi 0.4 a 0.9. Abych lépe viděl, co znamenají parametry v mých modelech, udělám transformaci tnox <- (nox-0.4)/0.1. Absolutní člen v mých modelech bude nyní udávat průměrnou cenu nemovitostí při koncentraci NOx = 0.4 · 10−7 (nikoli v 0) a parametr u nox bude udávat změnu ceny při nárůstu NOx o 0.1 · 10−7 (nikoli o 1). • První model, který vyzkouším, bude fit1 <- lm(medv~tnox). Dostanu 1 Michal
Kulich, KPMS MFF UK
Zobecněné lineární modely
Problém 1: Poznámky
Call: lm(formula = medv ~ tnox) Residuals: Min 1Q -13.691 -5.121
Median -2.161
3Q 2.959
Max 31.310
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 27.7795 0.6176 44.98 <2e-16 *** tnox -3.3916 0.3196 -10.61 <2e-16 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.323 on 504 degrees of freedom Multiple R-Squared: 0.1826, Adjusted R-squared: 0.181 F-statistic: 112.6 on 1 and 504 DF, p-value: < 2.2e-16 Vidím, že průměrná cena domu při koncentraci NOx = 0.4 jest $27780 a cena klesá o $3392 s nárůstem koncentrace o 0.1. Nicméně z výběrových kvantilů pro residua ve výšeuvedeném výpisu si hned všimnu, že residua jsou silně asymetrická. Totéž potvrdí obrázek qqnorm(resid(fit1)). • Zkusím ztransformovat cenu logaritmem. Dostanu model log(Y ) = β0 + β1 X + ε, čili Y = eβ0 eβ1 X eε . Mám tedy E ( Y | X = x ) = eβ0 eβ1 x E eε . Podělením E ( Y | X = x + 1 ) a E ( Y | X = x ) dostanu přesně eβ1 , takže 100(eβ1 − 1) můžu interpretovat jako percentuální přírůstek/úbytek E Y při změně X o jednu jednotku (tj. nárůstu koncentrace NOx o 0.1). Obecně nemohu tvrdit, že eβ0 je střední cena nemovitosti při koncentraci NOx = 0.4, ale za předpokladu symetrie rozdělení chyb ε lze říci, že eβ0 je mediánem ceny nemovitosti při koncentraci NOx = 0.4. Tak tedy fit2 <- lm(log(medv)~tnox): Call: lm(formula = log(medv) ~ tnox) Residuals: Min 1Q Median -1.17597 -0.19503 -0.03334
3Q 0.18223
Max 1.08159
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.31314 0.02610 126.92 <2e-16 *** tnox -0.18011 0.01351 -13.33 <2e-16 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3518 on 504 degrees of freedom Multiple R-Squared: 0.2607, Adjusted R-squared: 0.2592 F-statistic: 177.7 on 1 and 504 DF, p-value: < 2.2e-16 Se symetrií residuí jsme si dost pomohli, což potvrdí i qqnorm(fit2). Ještě ztransformuji parametry: > exp(coef(fit2)) (Intercept) tnox 27.4712390 0.8351754
2
Zobecněné lineární modely
Problém 1: Poznámky
30 20 10
medv
40
50
Obrázek 1: Porovnání lineárního a kubického modelu pro logaritmus ceny.
0
1
2
3
4
tnox
Odhadnutý medián ceny domu při koncentraci 0.4 je tedy $27471. Cena domu v průměru klesá o 16.5% při nárůstu koncentrace o 0.1. • Ještě se zabývejme nelineární transformací NOx . Zkusíme třeba polynom třetího řádu: > fit3 <- lm(log(medv)~poly(tnox,3)) > anova(fit2,fit3) Analysis of Variance Table Model 1: Model 2: Res.Df 1 504 2 502 ---
log(medv) ~ tnox log(medv) ~ poly(tnox, 3) RSS Df Sum of Sq F Pr(>F) 62.378 59.297 2 3.081 13.044 3.002e-06 ***
3
Zobecněné lineární modely Signif. codes:
Problém 1: Poznámky
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Jelikož původní model je vnořený do tohoto modelu, můžeme otestovat rozdíl mezi nimi F-testem. Ten je vysoce významný. Teď je nejlepší udělat si obrázek, abychom viděli, v čem se oba modely liší. Nakreslíme si predikce odezvy z obou modelů do jednoho obrázku (všimněte si použití funkce predict()): noxpts <- seq(min(tnox),max(tnox),length=300) newdata <- data.frame(tnox=noxpts) fitted2 <- exp(predict(fit2,newdata)) fitted3 <- exp(predict(fit3,newdata)) plot(tnox,medv) lines(noxpts,fitted2,lty=1) lines(noxpts,fitted3,lty=2) Výsledek je na obrázku 1. Vidíme, že (i) v rozmezí tnox od 0.5 do 3.5 je vztah v podstatě lineární a oba modely se liší jen málo; (ii) největší rozdíly mezi oběma jsou pro nejmenší a největší hodnoty znečištění; (iii) polynomiální model naznačuje, že pro znečištění menší než cca. 0.5 anebo větší než cca. 0.75 (na původní škále) je vliv koncentrace NOx na cenu minimální. Co se přesně na obou koncích děje, to je těžko soudit. Z obrázku 1 je vidět, že data o znečištění jsou podezřele seskupená — například hodnoty znečištění větší než 3 jsou prakticky diskrétní. To může být způsobeno metodikou měření (zaokrouhlování), nebo tím, že některé lokality jsou natolik blízko sebe, že jejich znečištění je prakticky stejné, nebo tím, že jeden měřicí přístroj udává hodnotu znečištění pro několik sousedních lokalit. V posledních dvou případech bychom měli problém s předpokladem nezávislosti — data ve skutečnosti přicházejí ve shlucích, ale o struktuře těchto shluků nemáme žádnou informaci. Závislost mezi pozorováními pak může způsobit zakřivení regresního vztahu při vysokých hodnotách znečištění. Budeme dál pokračovat v aplikaci lineárního modelu, ale měli bychom si uvědomovat, že předpoklady modelu nejspíš neplatí a být opatrní při interpretaci výsledků, které dostaneme. Nyní máme dvě možnosti: buď můžeme zvolit lineární závislost (a nelpět příliš na výsledcích pro oba extrémy NOx ), nebo přejít ke kubické závislosti (a mít potíže s vysvětlováním jeho parametrů). Je možné vymyslet i něco jiného (třeba spojitou po částech lineární křivku), ale vyberme si pro jednoduchost první variantu, lineární vztah. Vzhledem k tomu, že v modelu fit2 má znečištění vysoce významný vztah k ceně nemovitosti, můžeme si i při evidentním porušení předpokladu nezávislosti dovolit vyslovit dost jednoznačný závěr. • Odpovězme na otázku 1 takto: Cena nemovitosti statisticky významně souvisí se znečištěním ovzduší. Střední cena nemovitosti v různých lokalitách, které se liší v koncentraci NOx , klesá zhruba o 16.5% na každých 0.1 nárůstu koncentrace NOx . Vliv koncentrace NOx na velmi levné anebo velmi drahé nemovitosti může však být nižší než oněch průměrných 16.5%. • Otázka č. 2: Podívejme se nejprve, jak souvisí koncentrace NOx s ostatními veličinami: Toto jsou jejich korelace: crim zn indus chas nox rm age dis rad [1,] 0.421 -0.517 0.764 0.091 1 -0.302 0.731 -0.769 0.611 tax ptratio black lstat medv [1,] 0.668 0.189 0.348 0.591 -0.427 Vidíme, že vyšší znečištění může souviset s vyšší kriminalitou, nižším podílem velkých pozemků, vyšší industrialisací, vyšším stářím domů, větší blízkostí do centra, vyšší daňovou sazbou a vyšší mírou chudoby. Možná, že ve skutečnosti cenu pozemku ovlivňují jen tyto faktory, zatímco znečištění nehraje roli. Abychom zjistili, zdali tomu tak je, pokusíme se od vlivu znečištění odečíst vlivy těchto vedlejších matoucích (confounding) faktorů. V principu stačí sestavit model, který obsahuje kromě nox i ostatní potenciální vysvětlující veličiny, a podívat se, zdali i potom nox významně souvisí s cenou.
4
Zobecněné lineární modely
Problém 1: Poznámky
• Zkusíme model fitm1 <- lm(log(medv)~tnox+crim+zn+indus+chas+rm+age+dis+ rad+tax+ptratio+black+lstat) summary(fitm1) ... Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.9442422 0.1699369 23.210 < 2e-16 *** tnox -0.0784671 0.0152409 -5.148 3.80e-07 *** ... Člověk ihned vidí, že tento model vysvětluje cenu pozemku mnohem lépe a že řada z nově přidaných veličin má velmi úzký vztah k ceně pozemků. Těch pár, které nemají, můžeme ale nemusíme z modelu odstranit. Nás totiž nejvíc zajímá, co se stalo s koeficientem veličiny tnox. Ten je stále signifikantně různý od nuly, ale jeho hodnota se změnila. Hleďme: getci <- function(fit,var) { # get estimate and se for variable var a <- summary(fit)$coef[var,1:2] # get 95% confidence interv. ci <- rep(a[1],2)+c(-1,1)*a[2]*qt(0.975,fit$df.residual) names(ci) <- c("Lower","Upper") exp(c(a[1],ci)) } > getci(fit2,"tnox") Estimate Lower Upper 0.8351754 0.8132990 0.8576402 > getci(fitm1,"tnox") Estimate Lower Upper 0.9245324 0.8972575 0.9526365 Zavedli jsme novou funkci getci(fit,var), která vysaje z odhadnutého modelu fit výsledky pro veličinu ˆ var (zadat jako znakový řetězec v uvozovkách) a spočte eβ a 95% interval spolehlivosti pro eβ . Tuto funkci2 použijeme na fit2 a fitm1 a porovnáme výsledky pro tnox. Je vidět, že původní 16.5%-ní snížení průměrné ceny při vzrůstu koncentrace NOx o 0.1 (95%-ní interval spolehlivosti 14% – 19%; viz model fit2) se změnilo na 8.5%-ní snížení v modelu fitm1 (interval 5% – 10%). Je tedy vidět, že ostatní veličiny, které jsme vzali v úvahu, vysvětlují zhruba polovinu původně odhadnutéhu vlivu NOx , ale nikoli vliv celý. • Odpověď na otázku č. 2: Ostatní uvažované faktory vysvětlují zhruba polovinu vlivu koncentrace NOx na cenu nemovitosti. I když je vezmeme v úvahu, koncentrace NOx má stále negativní vztah k ceně nemovitosti a nárůst koncentrace o 0.1 vede ke snížení průměrné ceny o zhruba 8.5%. • Odpověď na otázku č. 3: Nelze vyloučit ani prokázat, že zlepšení čistoty ovzduší by vedlo ke zvýšení cen nemovitostí. I kdyby se tak stalo, jednalo by se pravděpodobně o zvýšení relativně malé. • Závěrečná poznámka: Kdybychom chtěli dělat věci pořádně, museli bychom řádně prozkoumat funkcionální vztah všech veličin v modelu fitm1 k odezvě log(medv) (tj., potenciální transformace regresorů) a museli bychom uvažovat možné interakce mezi nimi. 2 Jistě by šla použít i funkce confint v kombinaci s exp. getci je však snadno upravitelná i pro modely, na něž confint nefunguje.
5
