* Encoding: UTF-8.
* Budeme pracovat s daty za ČR.
USE ALL.
COMPUTE filter_$=(country=5).
VARIABLE LABELS filter_$ 'country=5 (FILTER)'.
VALUE LABELS filter_$ 0 'Not Selected' 1 'Selected'.
FORMATS filter_$ (f1.0).
FILTER BY filter_$.
EXECUTE.

*Máme 729 případů.
FREQUENCIES country.
 
* Naší závislou proměnnou bude čas užívání internetu timeuse.
EXAMINE timeuse
 /PLOT HISTOGRAM BOXPLOT.

DESCRIPTIVES VARIABLES=timeuse
  /STATISTICS=MEAN STDDEV MIN MAX.


* Naší nezávislou proměnnou bude socioekonomický status - ses.
FREQUENCIES ses.

*Jak se liší rozložení a průměry timeuse ve skupinách SES?.
MEANS TABLES=timeuse BY ses
  /CELLS=MEAN STDDEV COUNT.

EXAMINE VARIABLES=timeuse BY ses
  /PLOT BOXPLOT HISTOGRAM NPPLOT SPREADLEVEL(1)
  /COMPARE GROUPS
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Histogramy.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=timeuse ses MISSING=LISTWISE REPORTMISSING=NO
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
  SOURCE: s=userSource(id("graphdataset"))
  DATA: timeuse=col(source(s), name("timeuse"))
  DATA: ses=col(source(s), name("ses"), unit.category())
  GUIDE: axis(dim(1), label("Estimated minutes online each day"))
  GUIDE: axis(dim(2), label("Frequency Percent"))
  GUIDE: axis(dim(4), label("Socio-economic status - Household"), opposite())
  GUIDE: text.title(label("Simple Histogram of Estimated minutes online each day by ",
    "Socio-economic status - Household"))
  SCALE: cat(dim(4), include("1", "2", "3"))
  ELEMENT: interval(position(summary.percent.count(bin.rect(timeuse*1*1*ses, binCount(15)), 
    base.all())), shape.interior(shape.square))
END GPL.

*Průměr a 95% CI.
GRAPH
  /ERRORBAR(CI 95)=timeuse BY ses.

* Boxploty.
EXAMINE VARIABLES=timeuse BY ses
  /PLOT=BOXPLOT
  /STATISTICS=NONE
  /NOTOTAL.

*V obou chybí 10 lidí. Je to těch samých 10?. 
MVA VARIABLES=timeuse ses
  /MAXCAT=25
  /CATEGORICAL=ses
  /TTEST PROB PERCENT=0
  /TPATTERN PERCENT=0.
*Ne není. V ANOVA tak bude jen 709.

*Kdybychom nechtěli riskovat heteroskedascitu, můžřeme použít Welchovu korekci.
ONEWAY timeuse BY ses
  /STATISTICS DESCRIPTIVES HOMOGENEITY WELCH
  /PLOT MEANS
  /MISSING ANALYSIS.

*Můžeme reportovat, étu/omegu si musíme dopočítat:
 eta2 = BetweenSS / TotalSS
 omega2 = [BetweenSS - (DFbetween*MSwithin)] / (TotalSS + MSwithin)

* Existuje více "vzorečků", viz: 
http://www.jasonfinley.com/tools/OmegaSquaredQuickRef_JRF_3-31-13.pdf

* Kontrasty.
* Otestujme dvě hypotézy z přednášky.
    *H1: Děti z rodin s nízkým SES používají internet méně často než ostatní děti.
    *H2: Děti z rodin se středním SES používají internet méně často než děti z rodin s vysokým SES.
ONEWAY timeuse BY ses
  /STATISTICS DESCRIPTIVES HOMOGENEITY
  /CONTRAST -1 0.5 0.5 /*H1: Děti z rodin s nízkým SES používají internet méně často než ostatní děti
  /CONTRAST 0 -1 1  /*H2: Děti z rodin se středním SES používají internet méně často než děti z rodin s vysokým SES
  /PLOT MEANS
  /MISSING ANALYSIS.

* Kdybychom se zajímali o to, zda průměry skupin vykazují vykazují nějaký "trend"
  Polynomické konstrasty přímo neověřují, zda se průměry skupin liší, ale zda průměry skupin
  tvoří nějaký "trend", vzorec (např. přímku, parabolu).

ONEWAY timeuse BY ses
  /POLYNOMIAL=1 /* lineární kontrast
  /STATISTICS DESCRIPTIVES HOMOGENEITY
  /PLOT MEANS
  /MISSING ANALYSIS.

*Post-hoc testy?.
*Když už jsme použili kontrasty, není pro jejich použití sebemenší důvod.
*Takže jen abychom si je ukázali.
ONEWAY timeuse BY ses
  /STATISTICS DESCRIPTIVES   
  /POSTHOC = QREGW TUKEY
  /MISSING ANALYSIS.

*Porovnání s regresí.

*Dummy kódování.
    * .srovnání Low SES a High SES s referenční skupinou Middle SES. 
RECODE ses (1=1) (2=0) (3=0) INTO lowSES_dummy.
RECODE ses (1=0) (2=0) (3=1) INTO highSES_dummy.
EXECUTE.

REGRESSION
  /DESCRIPTIVES MEAN STDDEV CORR SIG N
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS CI(95) R ANOVA
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT timeuse
  /METHOD=ENTER lowSES_dummy highSES_dummy.

*Efektové kódování.
RECODE ses (1=1) (2= -1) (3=0) INTO lowSES_ec.
RECODE ses (1=0) (2=-1) (3=1) INTO highSES_ec.
EXECUTE.


REGRESSION
  /DESCRIPTIVES MEAN STDDEV CORR SIG N
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS CI(95) R ANOVA
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT timeuse
  /METHOD=ENTER lowSES_ec highSES_ec.

* Při efektovém kódování budou regresní koeficienty představovat rozdíl mezi průměrem dané skupiny (s kódem 1) 
  a celkovým neváženým průměrem (tj. průměr z průměrů jednotlivých skupin).
* Konstanta představuje celkový nevážený celkový průměr.

* Socio-economic status - Household	Mean
                   Mean	N
Low	148,84	58
Medium	128,97	407
High	115,70	244
Total_W	126,03	709
Total_UW    131.17        709
*
N_low/N_Medium = 0.142506
N_lhigh/N_Medium = 0.5995086.

*Vážené efektové kódování.
RECODE ses (1=1) (2= -0.142506) (3=0) INTO lowSES_wec.
RECODE ses (1=0) (2=-0.5995086) (3=1) INTO highSES_wec.
EXECUTE.

REGRESSION
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS R ANOVA
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT timeuse
  /METHOD=ENTER lowSES_wec highSES_wec.

* Při váženém efektovém kódování budou regresní koeficienty představovat rozdíl mezi průměrem dané skupiny (s kódem 1) 
  a celkovým (váženým) průměrem (průměr všech pozorování/případů bez ohledu na skupinu).
* Konstanta představuje tento celkový průměr.


*Přejděme k faktoriální anově.
*Nejprve si spočítáme one-way anovu pomocí příkazu UNIANOVA, abychom viděli, že výsledek je stejný.
UNIANOVA timeuse BY ses 
 /DESIGN = ses 
 /PRINT = DESCRIPTIVE 
 /PLOT PROFILE (ses) ERRORBAR = CI.

*Můžeme přidat další faktor - pohlaví.
*Máme ho, zdá se, pro všech 729 lidí, takže nám to neovlivní efektivní velikost vzorku.
FREQUENCIES child_gender.

MEANS TABLES=timeuse BY child_gender
  /CELLS=COUNT MEAN  STDDEV VAR.

*Dohromady i s SES.
CROSSTABS
  /TABLES= child_gender BY ses
  /FORMAT=AVALUE TABLES
  /CELLS=COUNT
  /COUNT ROUND CELL.

CTABLES
  /VLABELS VARIABLES=ses child_gender timeuse DISPLAY=LABEL
  /TABLE ses BY child_gender > timeuse [MEAN, STDDEV, COUNT]
  /CATEGORIES VARIABLES=ses child_gender ORDER=A KEY=VALUE EMPTY=INCLUDE
  /CRITERIA CILEVEL=95.

*Graficky. 
GRAPH
  /LINE(MULTIPLE)=MEAN(timeuse) BY ses BY child_gender
  /INTERVAL CI(95.0).

*Chci aspoň nakouknout na rozložení uvnitř skupin.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=ses timeuse child_gender MISSING=LISTWISE 
    REPORTMISSING=NO
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
  SOURCE: s=userSource(id("graphdataset"))
  DATA: ses=col(source(s), name("ses"), unit.category())
  DATA: timeuse=col(source(s), name("timeuse"))
  DATA: child_gender=col(source(s), name("child_gender"), unit.category())
  GUIDE: axis(dim(1), label("Socio-economic status - Household"))
  GUIDE: axis(dim(2), label("Estimated minutes online each day"))
  GUIDE: axis(dim(4), label("Gender of child?"), opposite())
  SCALE: cat(dim(1), include("1", "2", "3"))
  SCALE: linear(dim(2), include(0))
  SCALE: cat(dim(4), include("1", "2"))
  ELEMENT: point.dodge.symmetric(position(ses*timeuse*1*child_gender))
END GPL.


*To druhé BY nejde naklikat.
examine timeuse BY ses BY child_gender
 /PLOT HISTOGRAM BOXPLOT.

* A nyní si můžeme spočítat ANOVU.
UNIANOVA timeuse BY ses child_gender
 /DESIGN = ses child_gender ses*child_gender
 /PRINT = DESCRIPTIVE
  /PLOT PROFILE (ses*child_gender) ERRORBAR = CI.


UNIANOVA timeuse BY ses child_gender
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PLOT=PROFILE(ses*child_gender) TYPE=LINE ERRORBAR=CI MEANREFERENCE=NO YAXIS=0
  /PRINT ETASQ DESCRIPTIVE
  /CRITERIA=ALPHA(.05)
  /DESIGN=ses child_gender ses*child_gender.

*Protože máme v modelu více faktorů, můžeme požádat o EMM, 
tedy  odhady  toho, jaké by byly průměry ve skupinách podle jednoho faktoru,
kdyby v nich byly skupiny podle druhého faktoru zastoupeny stejně.
*Jsou to vlastně predikované skóry. 
UNIANOVA timeuse BY ses child_gender
 /DESIGN = ses child_gender ses*child_gender
 /PRINT = DESCRIPTIVE
  /PLOT PROFILE (ses*child_gender) ERRORBAR = CI
  /EMMEANS=TABLES(ses) 
  /EMMEANS=TABLES(child_gender). 

*Je dobré vědět, že defaultní model, který zahrnuje interakci, lze změnit.
UNIANOVA timeuse BY ses child_gender
 /DESIGN = ses child_gender 
 /PRINT = DESCRIPTIVE
  /PLOT PROFILE (ses*child_gender) ERRORBAR = CI
  /EMMEANS=TABLES(ses) 
  /EMMEANS=TABLES(child_gender). 

* A podívejme se na kontrasty. Způsob jejich zadávání se zde liší.
* Je dobré si uvědomit, že porovnáváme průměry podle jednoho faktoru korigované 
na možný vliv ostatních faktorů Estimated Marginal Means.   
UNIANOVA timeuse BY ses child_gender
  /CONTRAST(ses)=Helmert
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT ETASQ DESCRIPTIVE PARAMETER
  /CRITERIA=ALPHA(.05)  
  /DESIGN=ses child_gender ses*child_gender.

*Kontrast ručně.
UNIANOVA timeuse BY ses child_gender
  /CONTRAST(ses)=SPECIAL ( -1 .5 .5  0 -1  1)
  /CONTRAST(ses)=Helmert
  /EMMEANS=TABLES(ses) 
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT ETASQ DESCRIPTIVE PARAMETER
  /CRITERIA=ALPHA(.05)  
  /DESIGN=ses child_gender ses*child_gender.

*Příkaz UNIANOVA nabízí obrovské množství outputu. 
UNIANOVA timeuse BY ses child_gender
  /CONTRAST(ses)=Helmert
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PLOT=PROFILE(ses*child_gender) TYPE=LINE ERRORBAR=CI MEANREFERENCE=NO YAXIS=AUTO
  /EMMEANS=TABLES(ses) 
  /EMMEANS=TABLES(child_gender) 
  /PRINT MBP ETASQ DESCRIPTIVE PARAMETER HOMOGENEITY OPOWER
  /PLOT=SPREADLEVEL RESIDUALS
  /CRITERIA=ALPHA(.05)  
  /ROBUST=HC3     /*Žádost o analogii Welchova F
  /DESIGN=ses child_gender ses*child_gender.

*Nakonec ještě dvě jednotlivosti.
*Model s random faktorem.
USE ALL.
UNIANOVA timeuse BY country
 /DESIGN = country
 /RANDOM = country 
 /PRINT = DESCRIPTIVE
  /PLOT PROFILE (country) ERRORBAR = CI.

 UNIANOVA timeuse BY country
 /DESIGN = country
 /RANDOM = country 
 /PRINT = DESCRIPTIVE
  /EMMEANS=TABLES(country) COMPARE ADJ(SIDAK)
  /PLOT PROFILE (country) ERRORBAR = CI.

*ANCOVA.
UNIANOVA timeuse BY country WITH child_age
 /DESIGN = country child_age country*child_age
 /RANDOM = country 
 /PRINT = DESCRIPTIVE
  /PLOT PROFILE (country) ERRORBAR = CI.