clear, close all

%% Cvičení 11. Testy na vhodnost dat k regresi. Testy pro validaci regrese
% Table of Contents
% Příklad 1
% Příklad 2
% Příklad 3
% Příklad 4
% Samostatná práce v Graderu
%
% Jsou data vhodná k regresi?
% podle předpokladů zvolíme parametrický / neparametrický test na vhodnost k regresní analýze - viz. tabulka s testy.
% H0: veličiny jsou nezávislé
% HA:   nejsou nezávislé
% Pokud zamítáme H0 (p<\alpha ),  data jsou vhodná k lineární / nelineární regresi
% Regrese: Odhad regresních koeficientů, predikce, přímka / křivka
% Byla regrese vhodně zvolena?
% Validace lineární regrese - F-test:
% H0:  zvolená regrese nebyla správná
% HA:   byla správná
% Pokud zamítáme H0 (p<\alpha ),  regresní metoda byla správná
% Validace nelineární regrese - dwtest:
% H0:  zvolená regrese byla správná
% HA:   nebyla správná
% Pokud zamítáme H0 (p<\alpha ),  regresní metoda nebyla správná
%% Příklad 1
% Zajímá nás, zda u žáků základní školy je prospěch v matematice  závislý na době spánku. Údaje o bodovém hodnocení žáků z testu z matematiky a době spánku máme v tabulce. Testujte, zda jsou data  vhodná k regresi. V případě dat bez předpokladu normality použijte polynomiální regresi 3.řádu. Ověřte výsledky regrese.

spanek=[13.13058   10.906727   10.714454   10.43259   10.150883   9.7021626   8.2351257   8.1952563   7.9958005   7.8212473   7.6306907 6.9523974   6.8827789   6.6865477];
matematika=[99.972991   91.3185   90.76371   86.891455   86.626367   83.784329   82.982832   82.368872   81.914338   78.654401   75.383659 73.831872   68.342547   65.103527];

% Máme dva párové výběry
% Otestujeme data na vhodnost k regresi. Čím začneme?

%

% Jaký test na vhodnost dat k regresi použijeme?
%
% Oba výběry mají normalitu, chceme testovat data na vhodnost k regresi - použijeme Pearsonův test
% H0: doba spánku a výsledky z matematiky jsou  lineárně nezávislé
% HA:   nejsou  lineárně nezávislé

[rPear,p]=corr(spanek',matematika',"type","Pearson")

% Jsou data vhodná k lineární regresi?
% Provedeme lineární regresi $y=b_0+b_1x$

x=spanek; y=matematika;% x - nezávisle proměnná, y - závisle proměnná
mdl=fitlm(x,y,'linear')% % lineární regrese
b=mdl.Coefficients.Estimate % regresní koeficienty

% Pro vykreslení regresní přímky:

yp=predict(mdl,x');% data musí být ve sloupcích. Hodnoty na přímce - predikce
%yp5=predict(mdl,5) % predikce jenom pro 5 hodin spánku

figure
plot(x,y,'square')
hold on
plot(x,yp,'r')
xlabel('Doba spánku','Fontsize',12)
ylabel('Bodové hodnocení žáků','Fontsize',12)
legend('data','predikce','Location',"best")

%
% Pro validaci regrese použijeme F-test podílu vysvětleného a nevysvětleného rozptylu:
% H0:  zvolili jsme špatnou metodu
% HA:   zvolili jsme správnou metodu

% p-hodnota F-testu
p_F=coefTest(mdl)

% Byla regrese vhodně zvolena?
%% Příklad 2
% Zjišťujeme vztah mezi vývojem ceny benzínu natural 95 a řidičského kurzu skupiny B v období 2007-2017 v Kč. Otestujte, zda jsou  data vhodná k regresi a v případě dat bez předpokladu normality použijte polynomiální regresi 2.řádu. Následně ověřte správnost regrese.

benzin= [23.85  28.44  29.  30.  30.45   31.01   32.88   33.93   34.95  35.11   36.16 ];
kurz=[8658.17  9219.96  9146.95  9108.76  9206.27  9284.47  9316.35  9314.01  9316.69  9369.57   9416.43];

% Máme dva párové výběry
% Otestujeme data na vhodnost k regresi. Čím začneme?

%

% Jaký test na vhodnost dat k regresi použijeme?
%
% Jeden z výběrů nemá normalitu, chceme testovat data na vhodnost k regresi - použijeme Spearmanův test
% H0: cena benzínu a cena kurzu jsou  nezávislé
% HA:   nejsou nezávislé

[rSpear,p]=corr(benzin',kurz',"type","Spearman") % data musí být ve sloupcích

% Jsou data vhodná k regresi?
%
% Provedeme polynonomiální regresi $y=b_0+b_1x+b_2x^{2}$

x=benzin; y=kurz;% x - nezávisle proměnná, y - závisle proměnná
mdl2=fitlm(x,y,'poly2')% 2.řádu
b=mdl2.Coefficients.Estimate % regresní koeficienty

% Pro vykreslení regresní křivky:

yp2=predict(mdl2,x');% data musí být ve sloupcích. Hodnoty na křivce - predikce 
%yp18=predict(mdl2,18) % predikce pouze pro cenu 18Kč

figure
plot(x,y,'o')
hold on
plot(x,yp2,'r')
xlabel('Cena benzínu','Fontsize',12)
ylabel('Cena kurzu','Fontsize',12)
legend('data','predikce','Location',"best")

% Pro validaci regrese použijeme Durbinův-Watsonův test dwtest:
% H0:  zvolili jsme správnou metodu
% HA:   zvolili jsme špatnou metodu

p_dw=dwtest(mdl2)

% Byla regrese vhodně zvolena?
%
%% Příklad 3
% Hledáme, zda existuje souvislost mezi BMI, věkem a systolickým tlakem (závisle proměnná). Zeptali jsme 20 pacientů na jejich věk, BMI a změřili jim tlak. Na hladině významnosti 0,05 otestujte, zda jsou data vhodná k regresi. Pokud ano, proveďte regresi podle typu dat a poté ověřte správnost zvolené regrese.

tlak=[129 110 128 115 120 111 121 165 162 149 142 103 116 115 137 136 129 120 109 127];
BMI=[34.23 30.46 29.15 31.68 31.29 19.63 28.29 25.94 17.41 21.15 26.50 32.57 26.98 29.58 23.44 16.76 28.69 23.15 34.41 37.39];
vek=[36 19 50 57 44 21 59 45 50 52 49 42 36 25 63 61 64 41 26 63];

% Máme tři párové výběry - data se mohou hodit k vícenásobné lineární regresi
%                                                         $y=b_0+b_1x_{1}+b_2x_{2}$
% Otestujeme data na normalitu:

[h,p]=adtest(tlak)
[h,p]=adtest(BMI)
[h,p]=adtest(vek)

% Všechny výběry mají normalitu
% Nejdříve musíme testovat data na vhodnost k regresi - bereme  veličinu $y$ a každou veličinu $x$ zvlášť :

[rPear,p]=corr(tlak',BMI',"type","Pearson")
[rPear,p]=corr(tlak',vek',"type","Pearson")

% Jsou data vhodná k regresi?
% Naopak, všechny veličiny $x$ by měly být mezi sebou nezávislé:

[rPear,p]=corr(BMI',vek',"type","Pearson")

% Jsou  veličiny $x$ nezávislé?
%
% Provedeme vícenásobnou lineární regresi $y=b_0+b_1x_{1}+b_2x_{2}$

x1=BMI; x2=vek;y=tlak; % x1,x2 - nezávisle proměnné, y - závisle proměnná
mdl3=fitlm([x1; x2]',y)
yp2=predict(mdl3,[18 43])%predikce jenom pro BMI 18 a věk 43

% Byla regrese vhodně zvolena?
%% Příklad 4

T=readtable('ceny_CSU.xlsx',"Sheet","data_new","PreserveVariableNames",true) 

% Použijeme data ze statistického úřadu, vybereme si libovolné veličiny. Například:

[h,p]=adtest(T.("maso vepřové - bůček"))
[h,p]=adtest(T.("šunka vepřová"))

% Zamítáme normalitu druhého výběru - musíme použít neparametrický test:

[rSpear,p]=corr(T.("maso vepřové - bůček"),T.("šunka vepřová"),"type","Spearman")

% Jsou data vhodná k regresi? jaký druh regrese můžeme použít?

x=T.("maso vepřové - bůček"); y=T.("šunka vepřová");% x - nezávisle proměnná, y - závisle proměnná
mdl4=fitlm(x,y,'poly2')% 2.řádu
b=mdl4.Coefficients.Estimate % regresní koeficienty

% Jaký test validace regrese můžeme použít? Proč?

p_wd=dwtest(mdl4)

% Byla regrese vhodně zvolena?
%% Samostatná práce v Graderu
%
%
%