wyklad 6 z R 11

dane = read.csv2("wzrost_dziecka.csv")

summary(dane)

attach(dane)

(n=length(dane[,1]))

#zależnośc wzrost ojca/wzrost matki

#podejście macierzowe

X = cbind(rep(1, length(Ojciec)), Ojciec)

(XTX = t(X) %*% X)

(D=solve(XTX))

(a = D %*% t(X) %*% Matka)

#Błąd std. (odch std. reszt)

(blad = sqrt(sum((Matka - X %*% a)^2)/(n-2)))

(blad_ww = sqrt(D[1,1])*blad)

(blad_o = sqrt(D[2,2])*blad)

#czy reszty mają rozkład normalny?

shapiro.test(Matka - X %*% a)

#Test dla parametru przy wzroscie ojca (dla wyrazu wolnego analogicznie)

#H: a_1 = 0, K: a_1 <> 0

#stat. testowa

(t = a[2]/blad_o)

#p-wartość:

2*(1-pt(t, n-2))

#Wnioskowanie o regresji jako całości

(SST = var(Matka)*(n-1))

(SSR = sum((X %*% a - mean(Matka))^2))

(SSE = sum((Matka - X %*% a)^2))

#R^2

(R2 = SSR/SST)

#Test ANOVA

(StatF = (SSR/1)/(SSE/(n-2)))

#p-wartość

1-pf(StatF, 1, n-2)

#gotowe funkcje R:

plot(Ojciec, Matka)

reg1 = lm(Matka~Ojciec)

summary(reg1)

#korzysta m.in z:

anova(reg1)

#wykres linii regresji

ndane = data.frame(Ojciec = c(min(Ojciec), max(Ojciec)))

lines(ndane$Ojciec, predict(reg1, ndane), lwd = 2, col = "red")

#zależnośc wzrost matki/wzrost dziecka

plot(Matka, Dziecko)

reg2 = lm(Dziecko~Matka)

summary(reg2)

#wykres linii regresji

ndane = data.frame(Matka = c(min(Matka), max(Matka)))

lines(ndane$Matka, predict(reg2, ndane), lwd = 2, col = "purple")

#czy reszty mają rozkład normalny?

shapiro.test(reg2$residuals)

#Czyli wnioskowanie o regresji całkiem słuszne

#Wykres diagnostyczny

plot(reg2)

#Jeszcze jeden wykres + szacowanie p-ufności dla predykcji

ndane = data.frame(Matka = seq(min(Matka)-20, max(Matka)+20, length = 300))

pred = predict(reg2, ndane, interval = "prediction",level =0.95)

plot(Matka, Dziecko)

lines(ndane$Matka, pred[,1], lwd = 2, col = "purple")

lines(ndane$Matka, pred[,2], lwd = 1, col = "red")

lines(ndane$Matka, pred[,3], lwd = 1, col = "red")

plot(ndane$Matka, pred[,1], lwd = 2, col = "purple", type = "l")

points(Matka, Dziecko)

lines(ndane$Matka, pred[,2], lwd = 1, col = "red")

lines(ndane$Matka, pred[,3], lwd = 1, col = "red")

#niektóre inne modele regresji daje się przez przekształcenia doprowadzić do liniowego

#ale trzeba ostrożniej podchodzić do oceny parametrów

reg3 = lm(Dziecko~I(Matka^2)+Matka)

summary(reg3)

#wykres linii regresji

ndane = data.frame(Matka = seq(min(Matka), max(Matka), length = 300))

lines(ndane$Matka, predict(reg3, ndane), lwd = 2, col = "violet")

#nieliniowe modele regresji - estymacja numeryczna

reg4 = nls(Dziecko~a*Matka^2+b*Matka+d, start=list(a=0, b=0, d=0))

summary(reg4)

#wykres linii regresji - oczywiście wyszło to samo

lines(ndane$Matka, predict(reg3, ndane), lwd = 2, col = "pink")

#regresja od większej ilości zmiennych (choć coś już się pojawiło)

#zależnośc wzrost dziecka/wzrost rodziców

reg5 = lm(Dziecko~Matka+Ojciec)

summary(reg5)

#super model!

#wykres - tu będzie trudniej - pseudo 3d

#matka od 155 do 180, ojciec od 165 do 195

mv = seq(155, 180, by = 5)

ov = seq(165, 195, by = 5)

#funkcja regresji

fun = function(x, y) reg5$coef[1] + x * reg5$coef[2] + y * reg5$coef[3]

dv = outer(mv, ov, fun)

wykres = persp(x=mv, y = ov, z = dv, xlab = "Matka", ylab = "Ojciec",

zlab = "Dziecko", ticktype = "detailed")

#nieładny punkt widzenia, to może to:

wykres = persp(x=mv, y = ov, z = dv, xlab = "Matka", ylab = "Ojciec",

zlab = "Dziecko", ticktype = "detailed", theta = 35, phi = 25,

cex.axis = 0.7)

#dane empiryczne:

points(trans3d(Matka, Ojciec, Dziecko, wykres), col="red", lwd = 2, pch = 20)

#to jeszcze kreski łączące wartości empiryczne i teoretyczne, żeby

#rozpoznać gdzie faktycznie leżą te punkty

segments(trans3d(Matka, Ojciec, Dziecko, wykres)$x,

trans3d(Matka, Ojciec, Dziecko, wykres)$y,

trans3d(Matka, Ojciec, predict(reg5), wykres)$x,

trans3d(Matka, Ojciec, predict(reg5), wykres)$y)

#to jeszcze dwa interesujące modele:

# Dziecko = a*(Matka+Ojciec)+b

reg6 = lm(Dziecko~I(Matka+Ojciec))

summary(reg6)

#ten jest równie dobry!

#wykres 2d!

plot((Matka+Ojciec)/2, Dziecko)

ndane = data.frame(Ojciec = c(min(Ojciec), max(Ojciec)), Matka = c(min(Matka), max(Matka)))

lines((ndane$Ojciec+ndane$Matka)/2, predict(reg6, ndane), lwd = 2, col = "green")

#Diagnostyka

plot(reg6)

#wykres z p.ufnosci

ndane = data.frame(Ojciec = seq(min(Ojciec)-10, max(Ojciec)+10, length = 300),

Matka = seq(min(Matka)-10, max(Matka)+10, length = 300))

pred = predict(reg6, ndane, interval = "prediction",level =0.95)

plot((ndane$Ojciec+ndane$Matka)/2, pred[,1], type = "l", lwd = 2, col = "green")

points((Matka+Ojciec)/2, Dziecko)

lines((ndane$Ojciec+ndane$Matka)/2, pred[,2], lwd = 1, col = "red")

lines((ndane$Ojciec+ndane$Matka)/2, pred[,3], lwd = 1, col = "red"

#model bez wyrazu wolnego:

reg7 = lm(Dziecko~Matka+Ojciec+0)

summary(reg7)

#uwaga na R^2 - nie ma sensu

detach(dane)

#ciąg (jeszcze) dalszy nastąpi

…..

#regresja logistyczna - dane Raty.csv

raty = read.csv2("Raty.csv")

attach(raty)

#dla wygody (wykresy) przekształćmy zmienną Raty na wartości 0 i 1.

(raty01 = c(Raty) - 1)

#regresja logistyczna - od wydatku

regl1 = glm(raty01 ~ Wydatek, family = "binomial")

summary(regl1)

plot(Wydatek, raty01)

#wykres linii regresji

x = seq(min(Wydatek), max(Wydatek), length = 300)

ndane = data.frame(Wydatek = x)

lines(x, 1/(1+exp(-predict(regl1, ndane))), lwd = 2, col = "red")

#alternatywnie

y = 1/(1 + exp(-regl1$coef[1] - regl1$coef[2] * x))

lines(x, y, col = "black")

#można zaobserwować zależność: im większy wydatek tym (trochę) częstsza

#decyzja o zakupie na raty.

#regresja logistyczna od dwóch zmiennych

regl2 = glm(raty01 ~ Wydatek+Dochody, family = "binomial")

summary(regl2)

#wykres powierzchni regresji - pseudo 3d

#wydatek od 200 do 2700, dochody od 1000 do 4100

wydv = seq(200, 2700, by = 100)

docv = seq(1000, 4100, by = 100)

#funkcja regresji

f_l = function(x, y) 1/(1+exp(-regl2$coef[1] - x * regl2$coef[2] - y * regl2$coef[3]))

ratv = outer(wydv, docv, f_l)

wykres = persp(x=wydv, y = docv, z = ratv, xlab = "Wydatek", ylab = "Dochody",

zlab = "Raty", ticktype = "detailed", theta = 145, phi = 25,

cex.axis = 0.7, col = "yellow", shade = 0.6)

#dane empiryczne:

points(trans3d(Wydatek, Dochody, raty01, wykres), col="brown", lwd = 2, pch = 20)

#to jeszcze kreski łączące wartości empiryczne i teoretyczne, żeby

#rozpoznać gdzie faktycznie leżą te punkty

segments(trans3d(Wydatek, Dochody, raty01, wykres)$x,

trans3d(Wydatek, Dochody, raty01, wykres)$y,

trans3d(Wydatek, Dochody, f_l(Wydatek, Dochody), wykres)$x,

trans3d(Wydatek, Dochody, f_l(Wydatek, Dochody), wykres)$y)

#a może lepiej rzut na wykres dwuwymiarowy?

logit = predict(regl2)

plot(logit, raty01)

x = seq(min(logit), max(logit), length = 300)

y = 1/(1+exp(-x))

lines(x, y, lwd = 2, col = "blue2")

#widać, że model jest niezły

detach(raty)

#nieparametryczne metody regresji - z powrotem dane wzrost_dziecka.csv

attach(dane)

#regresja lokalnie wielomianowa

reg8 = loess(Dziecko~Matka)

summary(reg8)

plot(Matka, Dziecko)

ndane = data.frame(Matka = seq(min(Matka), max(Matka), length = 300))

lines(ndane$Matka, predict(reg8, ndane), lwd = 2, col = "green")

lines(ndane$Matka, predict(loess(Dziecko~Matka, span = 0.6, degree = 1), ndane), lwd = 2, col = "blue")

lines(ndane$Matka, predict(loess(Dziecko~Matka, span = 0.4, degree = 2), ndane), lwd = 2, col = "red")

#gładkie funkcje sklejane

#zależnośc Dziecko~Matka

#domyślne parametry

reg9 = smooth.spline(Matka, Dziecko)

print(reg9)

#wykres

plot(Matka, Dziecko)

lines(predict(reg9, x = seq(min(Matka), max(Matka), length = 400)), lwd = 2, col = "red")

lines(ndane$Matka, predict(reg8, ndane), lwd = 2, col = "green")

#odp. R^2

1-sum((Dziecko-predict(reg9, x = Matka)$y)^2)/sum((Dziecko-mean(Dziecko))^2)

#to jeszcze dwie linie (parametr spar)

#spar = 0 => lambda (prawie) równa 0

reg10 = smooth.spline(Matka, Dziecko, spar = 0)

lines(predict(reg10, x = seq(min(Matka), max(Matka), length = 400)), lwd = 2, col = "blue")

print(reg10)

#spar 1 => lamda duża

reg11 = smooth.spline(Matka, Dziecko, spar = 1)

lines(predict(reg11, x = seq(min(Matka), max(Matka), length = 400)), lwd = 2, col = "brown")

#praktycznie to samo, co regresja liniowa

print(reg11

detach(dane)


Wyszukiwarka

Podobne podstrony:
wyklad 11
WYKŁAD 11 SPS 2 regulatory 0
wyklad 11 toksyczno niemetali
BUD OG wykład 11 3 Geosyntetyki
Psychometria 2009, Wykład 11, Inwentarz MMPI
BUD OG wykład 11 1 Tworzywa sztuczne
Wyklad 11 2010
Wyklad 2 11
F II wyklad 11 30 04 12
chem wykład 11
Chemia fizyczna wykład 11
6 Miedzynarodowy transfer wyklad 11 04 2012 id 43355
Socjologia - wykład 11, geografia UJ, socjologia, wykłady 2010
Wykład 11.01.15 - Audiologia, Logopedia - podyplomowe, I sem - Audiologia
005 Historia sztuki wczesnochrześcijańskiej i bizantyjskiej, wykład, 11 09

więcej podobnych podstron