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)