ZADANIE 9

Kristína Mečiarová

---

V tomto zadaní budeme opäť pracovať s dátovým súborom uscrime z knižnice SMS data. Súbor obsahuje 50 pozorovaní a 11 premenných s informáciami o počte výskytov rôznych trestných činov v rámci konkrétnych regiónov USA, ich rozlohe a populácii.

1. LDA klasifikácia

Obrázok nižšie zobrazuje, ako bola klasifikované odozva. Na základe histogramov by sme priradili pozorovanie do západného regiónu, ak by skalárny súčin vektoru parametrov a vektoru koeficientov z lineárnej diskriminačnej analýzy nadobúdal skôr záporné hodnoty a pre kladné hodnoty tohto skalárneho súčinu by sme skôr usudzovali, že sa jedná o východný región.

library("Epi") #for relevel
data(uscrime, package = "SMSdata")
data<-transform(uscrime,newreg = Relevel(reg, list("West" = c("West","South"),"East" =c("Northeast","Midwest"))))

library(MASS)
lda1 <- lda(newreg ~ land.area + popu.1985 + murder + rape+ robbery + assault + burglary + larceny + autotheft, data = data)

pred <- predict(lda1)
par(mar=c(4,4,4,4))
ldahist(data = pred$x[,1], g=data$newreg)

Keďže máme 2 kategórie, máme iba 1 smer maximalizujúci variabilitu medzi skupinamu a variabilitu vrámci skupín. Nasledujúci obrázok zobrazuje, ako sú skupiny oddelené vzhľadom k tomuto smeru.

tData <- cbind(as.matrix(data[, c("land.area", "popu.1985", "murder", "rape", "robbery","assault","burglary","larceny","autotheft")]) %*% lda1$scaling, data$newreg)
m1 <- mean(tData[tData[,2] == 1, 1])
m2 <- mean(tData[tData[,2] == 2, 1])
plot(tData[,1], pch = 22, bg = c("yellow", "green")[unclass(as.factor(data$newreg))], ylab = " LDA direction",xlab="States")
abline(h=m1, lw=2, col = "yellow")
abline(h=m2, lw=2, col = "green")
legend("topright",legend=levels(data$newreg), pch=21,col=c("yellow", "green"),cex=0.8)

Pozrieme sa aj na presnosť predikcie. Nasledujúca tabuľka zobrazuje maticu zámien pre LDA a na obrázku sú vykreslené jednotlivé priradené pozorovania. Ako vidíme aj v tabuľke, 2 pozorovania z regiónu West boli chybne priradené do regiónu East a 2 pozorovania z regiónu East boli chybne priradené do regiónu West (farebné odlíšenie). Presnosť predikcie pre jednotlivé regióny je vypísaná pod tabuľkou.

(conf<-table(predict(lda1,type="class")$class,data$newreg))

##       
##        West East
##   West   27    2
##   East    2   19

diag(conf) / rowSums(conf)

##      West      East 
## 0.9310345 0.9047619

plot(pred$x[,1], pred$class, col=as.numeric(data$newreg)+10,ylab="Region",xlab="Predicted values")
legend("bottomright",legend=levels(data$newreg), pch=21,col=unique(as.numeric(data$newreg))+10,cex=0.8)

2. QDA klasifikácia

Ako alternatívnu analýzu použijeme QDA, tentokrát však iba pre premenné population a rape.

qda1<-qda(newreg ~rape +popu.1985, data = data)
lda2<-lda(newreg ~rape +popu.1985, data = data)

variab.gr = expand.grid(seq(min(data$popu.1985),max(data$popu.1985),len=1000),seq(min(data$rape),max(data$rape),len=1000))
names(variab.gr) = c("popu.1985","rape")

par(mfrow=c(1,2))
cols <- c("yellow","green")
variab.pop = unique(variab.gr$popu.1985)
variab.rap = unique(variab.gr$rape)

plot(data$popu.1985,data$rape,pch=20 + as.numeric(data$newreg), bg = cols[data$newreg], main = "Linear Discriminant Analysis",xlab="Population",ylab="Rape")
contour(variab.pop,variab.rap,
        matrix(as.numeric(predict(lda2,new=variab.gr)$class),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)

plot(data$popu.1985,data$rape,pch=20 + as.numeric(data$newreg), bg = cols[data$newreg], xlab="Population",ylab="Rape",main = "Quadratic Discriminant Analysis")
contour(variab.pop,variab.rap,
        matrix(as.numeric(predict(qda1,new=variab.gr)$class),1000,1000),
        nlev=1,drawlabels=FALSE,add=T)

par(mfrow=c(1,1))

Prvá tabuľka zachytáva predikcie pre LDA a druhá tabuľka pre QDA. Ako vidíme, tieto výsledky sú veľmi podobné, čo môžeme pozorovať aj na obrázkoch nižšie.

(conf2<-table(predict(lda2,type="class")$class,data$newreg))

##       
##        West East
##   West   24    7
##   East    5   14

(conf1<-table(predict(qda1,type="class")$class,data$newreg))

##       
##        West East
##   West   23    6
##   East    6   15

pred1 <- predict(qda1)
pred2 <- predict(lda2)

par(mfrow=c(1,2))
plot(pred2$x[,1], pred2$class, col=as.numeric(data$newreg)+10,main="LDA",ylab="Region",xlab="Predicted values")
legend("topleft",legend=levels(data$newreg), pch=21,col=unique(as.numeric(data$newreg))+10,cex=0.5)

plot(pred1$posterior[,2], pred1$class, col=as.numeric(data$newreg)+10,main="QDA",ylab="Region",xlab="Predicted values")
legend("topleft",legend=levels(data$newreg), pch=21,col=unique(as.numeric(data$newreg))+10,cex=0.5)

par(mfrow=c(1,1))