U 15. delu smo pokazali kako oceniti jednostavni uopšteni linearni model na binarnim podacima pomoću glm() naredbe. Nastavljamo sa istom glm() naredbom (modeliranje promenljive na težinu i radni takt motora) na istom skupu podataka mtcars koji je ugrađen u R.
model <- glm(formula = vs ~ wt + disp, data = mtcars, family = binomial)
summary(model)

Vidimo da uticaj težine pozitivno utiče na vs, dok radni takt ima mali negativan efekat. Takođe vidimo da koeficijent težine nije statistički značajan (p-vrednost < 0.05) dok je koeficijent radnog takta značajan. Kasnije ćemo videti načine poboljšanja našeg modela.
Međutim, ocene (koeficijenti objašnjavajućih promenljivih težine i radnog takta) sada se nalaze u jedinicama zvanim logit. Definisaćemo logit u kasnijem blogu.
Vidimo reč odstupanje (eng. Deviance) dva puta u modelu. Odstupanje je mera prilagođenosti modela podacima. R izveštava o dva oblika odstupanja – nulto odstupanje i rezidualno odstupanje. Nulto odstupanje pokazuje koliko model koji uključuje samo slobodni član, tj. sveopštu srednju vrednost (engl. intercept, i.e.grand mean) dobro predviđa promenljivu odziva. U našem primeru imamo vrednost hi-kvadrat statistike 43.9 sa 31 stepenom slobode. Ova vrednost ukazuje na malu prilagođenost, tj. značajnu razlika između ocenjenih i posmatranih vrednosti. Uključujući nezavisne promenljive (težine i radnog takta), smanjuje odstupanje na 21.4 sa 29 stepeni slobode, tako da ova vrednost hi-kvadrat statistike pokazuje značajno smanjenje odstupanja. Rezidualno odstupanje je smanjeno na 21.4 sa gubitkom od tri stepena slobode. Imamo veoma malu i samim tim i veoma statistički značajnu p-vrednost za model.
1 - pchisq(21.4, 3)

Primenimo Hosmer i Lemeshow (GOF) test testirajući hipotezu da se model dobro prilagođava podacima.
library(ResourceSelection)
hoslem.test(mtcars$vs, fitted(model))

Na žalost, moramo odbaciti nultu hipotezu, tj. naš model, stoga što prilagođene vrednosti su značajno različite od odgovarajućih podataka.
plot(sort(model$fitted), pch = 16, cex.lab = 1.3, ylim = c(0, 1))
points(sort(mtcars$vs), pch = 19, col = "red", cex = 1.3)

Kreirajmo drugi grafik, spajajući ovoga puta prilagođene vrednosti pravim linijama.
plot(sort(mtcars$vs), pch = 19, col = "red", cex = 1.3, ylab = "Proporcija")
points(sort(model$fitted), ty = "o", pch = 16, cex.lab = 1.3, ylim = c(0, 1))

Naredba glm() je generisala neku vrstu funkcije koja povezuje observacije vs = 0 sa vs = 1.
Možemo sada nacrtati grafik za svaku objašnjavajuću promenljivu. Prvo ćemo prilagoditi glm() promenljivoj težina.
model_težina <- glm(formula = vs ~ wt, data = mtcars, family = binomial)
summary(model_težina)

range(mtcars$wt)

xtežina <- seq(0, 6, 0.01)
ytežina <- predict(model_težina, list(wt = xtežina), type = "response")
plot(mtcars$wt, mtcars$vs, pch = 16, xlab = "TEŽINA (g)", ylab = "VS")
lines(xtežina, ytežina)

model_disp <- glm(formula= vs ~ disp, data = mtcars, family = binomial)
summary(model_disp)

range(mtcars$disp)

xdisp <- seq(0, 500, 0.01)
ydisp <- predict(model_disp, list(disp = xdisp), type = "response")
plot(mtcars$disp, mtcars$vs, pch = 16, xlab = "RADNI TAKT (kubni inči)", ylab = "VS")
lines(xdisp, ydisp)

To nije bilo tako teško! U 17. delu ću dati više objašnjenja rezultata koji smo dobili funkcijom glm().
Vidimo se kasnije!
David

Dodatak: Korišćeni R kodovi

# Ocena Uopštenog Linearnog Modela i prikaz rezultata.
model <- glm(formula = vs ~ wt + disp, data = mtcars, family = binomial) 
summary(model)

# Izračunavanje hi-kvadrat statistike.
1 - pchisq(21.4, 3)

# Učitavanje ResourceSelection paketa. 
library(ResourceSelection)

# Primena Hosmer i Lemeshow testa.
hoslem.test(mtcars$vs, fitted(model))

# Grafik ocenjenog modela. 
plot(sort(model$fitted), pch = 16, cex.lab = 1.3, ylim = c(0, 1))
points(sort(mtcars$vs), pch = 19, col = "red", cex = 1.3)

# Alternativni grafik ocenjenog modela. 
plot(sort(mtcars$vs), pch = 19, col = "red", cex = 1.3, ylab = "Proporcija")
points(sort(model$fitted), ty = "o", pch = 16, cex.lab = 1.3, ylim = c(0, 1))

# Ocena Uopštenog Linearnog Modela i prikaz rezultata.
model_težina <- glm(formula = vs ~ wt, data = mtcars, family = binomial)
summary(model_težina)

# Calculate range of values. 
range(mtcars$wt)

# Grafik ocenjenog modela. 
xtežina <- seq(0, 6, 0.01)
ytežina <- predict(model_težina, list(wt = xtežina), type = "response")
plot(mtcars$wt, mtcars$vs, pch = 16, xlab = "TEŽINA (g)", ylab = "VS")
lines(xtežina, ytežina)

# Ocena Uopštenog Linearnog Modela i prikaz rezultata.
model_disp <- glm(formula = vs ~ disp, data = mtcars, family = binomial)
summary(model_disp)

# Grafik ocenjenog modela. 
xdisp <- seq(0, 500, 0.01)
ydisp <- predict(model_disp, list(disp = xdisp), type = "response")
plot(mtcars$disp, mtcars$vs, pch = 16, xlab = "RADNI TAKT (kubni inči)", ylab = "VS")
lines(xdisp, ydisp)