% Options for packages loaded elsewhere
\PassOptionsToPackage{unicode}{hyperref}
\PassOptionsToPackage{hyphens}{url}
%
\documentclass[
]{ctexart}
\usepackage{amsmath,amssymb}
\usepackage{iftex}
\ifPDFTeX
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
  \usepackage{textcomp} % provide euro and other symbols
\else % if luatex or xetex
  \usepackage{unicode-math} % this also loads fontspec
  \defaultfontfeatures{Scale=MatchLowercase}
  \defaultfontfeatures[\rmfamily]{Ligatures=TeX,Scale=1}
\fi
\usepackage{lmodern}
\ifPDFTeX\else
  % xetex/luatex font selection
\fi
% Use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\IfFileExists{microtype.sty}{% use microtype if available
  \usepackage[]{microtype}
  \UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\makeatletter
\@ifundefined{KOMAClassName}{% if non-KOMA class
  \IfFileExists{parskip.sty}{%
    \usepackage{parskip}
  }{% else
    \setlength{\parindent}{0pt}
    \setlength{\parskip}{6pt plus 2pt minus 1pt}}
}{% if KOMA class
  \KOMAoptions{parskip=half}}
\makeatother
\usepackage{xcolor}
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\usepackage{graphicx}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
% Set default figure placement to htbp
\makeatletter
\def\fps@figure{htbp}
\makeatother
\setlength{\emergencystretch}{3em} % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{5}
\usepackage[most]{tcolorbox}
\definecolor{boxidea}{HTML}{6c8aa5}
\definecolor{boxwarn}{HTML}{b58266}
\definecolor{boxtip}{HTML}{8aa085}
\newtcolorbox{concept}[1][]{enhanced, breakable, colback=boxidea!7, colframe=boxidea, coltitle=white, fonttitle=\bfseries, title=#1, arc=2pt, left=8pt, right=8pt, top=4pt, bottom=4pt}
\newtcolorbox{pitfall}[1][]{enhanced, breakable, colback=boxwarn!7, colframe=boxwarn, coltitle=white, fonttitle=\bfseries, title=#1, arc=2pt, left=8pt, right=8pt, top=4pt, bottom=4pt}
\newtcolorbox{keypoint}[1][]{enhanced, breakable, colback=boxtip!7, colframe=boxtip, coltitle=white, fonttitle=\bfseries, title=#1, arc=2pt, left=8pt, right=8pt, top=4pt, bottom=4pt}
\ifLuaTeX
  \usepackage{selnolig}  % disable illegal ligatures
\fi
\IfFileExists{bookmark.sty}{\usepackage{bookmark}}{\usepackage{hyperref}}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\urlstyle{same}
\hypersetup{
  pdftitle={Nested ANOVA and MANOVA},
  pdfauthor={Rui Zhou},
  hidelinks,
  pdfcreator={LaTeX via pandoc}}

\title{Nested ANOVA and MANOVA}
\author{Rui Zhou}
\date{}

\begin{document}
\maketitle

{
\setcounter{tocdepth}{2}
\tableofcontents
}
\hypertarget{ux8fd9ux7ae0ux5728ux8bb2ux4ec0ux4e48}{%
\section{这章在讲什么}\label{ux8fd9ux7ae0ux5728ux8bb2ux4ec0ux4e48}}

普通单因素 ANOVA
假设组与组之间独立、组内观测也独立。但心理学实验里数据经常\textbf{有层级结构}：每个学校抽几个班，每个班抽几个学生------「学生」嵌套在「班」里，「班」嵌套在「学校」里。如果不显式建模这种嵌套，就会低估误差、把
\(p\) 值算小、得出虚高的显著性。这一章讲两件事：\textbf{嵌套
ANOVA}（一个因变量，但变量有嵌套层级）和
\textbf{MANOVA}（多个因变量同时考察组间差异）。

\begin{concept}[两个核心区别]
\textbf{Nested ANOVA}：误差项要拆——「组间」的检验要用「子组间」的方差作分母，而不是「子组内」。R 里 \texttt{aov(DV \~{} group + Error(group/subgroup))} 这种 \texttt{Error()} 项就是干这件事。

\textbf{MANOVA}：把多个 DV 看成一个向量 $\mathbf{Y}_i$，检验「组别」是否对 $\mathbf{Y}$ 的多维分布有影响。统计量从 $F$ 换成 Pillai's trace、Wilks' $\Lambda$、Hotelling-Lawley、Roy's largest root——这四个统计量各有侧重，Pillai 最稳健（对协方差齐性违反不敏感），是默认选择。MANOVA 显著之后再做 univariate ANOVA 看是哪个 DV 在驱动差异。
\end{concept}

\hypertarget{nested-anova}{%
\section{Nested ANOVA}\label{nested-anova}}

\hypertarget{section}{%
\subsection{0}\label{section}}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(tidyverse)}
\FunctionTok{library}\NormalTok{(car)}

\NormalTok{nestedData }\OtherTok{=}\NormalTok{ readxl}\SpecialCharTok{::}\FunctionTok{read\_excel}\NormalTok{(}\StringTok{"nested anova.xlsx"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\hypertarget{visualization}{%
\subsection{visualization}\label{visualization}}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# we should eyeball the data first to make sense of the test.}
\FunctionTok{library}\NormalTok{(ggplot2)}
\CommentTok{\# create boxplots to visualize the data}
\FunctionTok{ggplot}\NormalTok{(nestedData, }\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\FunctionTok{factor}\NormalTok{(location), }\AttributeTok{y=}\NormalTok{DV)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_boxplot}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{r-anova-manova_files/figure-latex/unnamed-chunk-2-1.pdf}

\hypertarget{nested-anova-1}{%
\subsection{nested ANOVA}\label{nested-anova-1}}

Theoretically, the order of computing the variance sum of squares
matters, thus appointing the
\texttt{summation\ =\ \textquotesingle{}I\textquotesingle{}}.

However, here the summation type does not matter.

Note that the result of the main effect of the group is wrong. We can
manually make it straight from the SS of the group and the subgroup.

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# note that nested anova can use error type I and III.}
\CommentTok{\# aov() has a default type I to process factors in an sequential order}
\CommentTok{\# here type I and III don\textquotesingle{}t differ, since effectively we have no orders}
\NormalTok{aov\_model }\OtherTok{=} \FunctionTok{aov}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(method)}\SpecialCharTok{/}\FunctionTok{factor}\NormalTok{(location), }
                \AttributeTok{data =}\NormalTok{ nestedData,}\AttributeTok{summation=}\StringTok{"I"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'summation' will be disregarded
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary}\NormalTok{(aov\_model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                 Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(method)                   2  665.7   332.8   255.7 1.45e-10 ***
## factor(method):factor(location)  9 1720.7   191.2   146.9 6.98e-11 ***
## Residuals                       12   15.6     1.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{aov\_model }\OtherTok{=} \FunctionTok{aov}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(method)}\SpecialCharTok{/}\FunctionTok{factor}\NormalTok{(location), }
                \AttributeTok{data =}\NormalTok{ nestedData,}\AttributeTok{summation=}\StringTok{"III"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'summation' will be disregarded
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary}\NormalTok{(aov\_model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                 Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(method)                   2  665.7   332.8   255.7 1.45e-10 ***
## factor(method):factor(location)  9 1720.7   191.2   146.9 6.98e-11 ***
## Residuals                       12   15.6     1.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\hypertarget{error-adjustment}{%
\subsubsection{error adjustment}\label{error-adjustment}}

An alternatDVe refinement is to adjust the error term.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{aov.result2 }\OtherTok{=} \FunctionTok{aov}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(method) }\SpecialCharTok{+} 
                    \FunctionTok{Error}\NormalTok{(}\FunctionTok{factor}\NormalTok{(method)}\SpecialCharTok{:}\FunctionTok{factor}\NormalTok{(location)),}
                  \AttributeTok{data =}\NormalTok{ nestedData)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in aov(DV ~ factor(method) + Error(factor(method):factor(location)), :
## Error() model is singular
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary}\NormalTok{(aov.result2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Error: factor(method):factor(location)
##                Df Sum Sq Mean Sq F value Pr(>F)
## factor(method)  2  665.7   332.8   1.741   0.23
## Residuals       9 1720.7   191.2               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 12  15.62   1.302
\end{verbatim}

\hypertarget{comparison-with-one-way}{%
\subsubsection{comparison with one-way}\label{comparison-with-one-way}}

If we wrongly ignore the nested factor and use one-way ANOVA, the
results would be different

Note that there's only one group factor, therefore, the type of
summation does not matter.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{aov\_model }\OtherTok{=} \FunctionTok{aov}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(method), }
                \AttributeTok{data =}\NormalTok{ nestedData,}\AttributeTok{summation=}\StringTok{"III"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'summation' will be disregarded
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary}\NormalTok{(aov\_model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                Df Sum Sq Mean Sq F value Pr(>F)  
## factor(method)  2  665.7   332.8   4.026 0.0331 *
## Residuals      21 1736.3    82.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{aov\_model }\OtherTok{=} \FunctionTok{aov}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(method), }
                \AttributeTok{data =}\NormalTok{ nestedData)}
\FunctionTok{summary}\NormalTok{(aov\_model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                Df Sum Sq Mean Sq F value Pr(>F)  
## factor(method)  2  665.7   332.8   4.026 0.0331 *
## Residuals      21 1736.3    82.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\hypertarget{manova}{%
\section{MANOVA}\label{manova}}

\hypertarget{section-1}{%
\subsection{0}\label{section-1}}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(}\StringTok{"readxl"}\NormalTok{)}
\FunctionTok{library}\NormalTok{(}\StringTok{"rstatix"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Attaching package: 'rstatix'
\end{verbatim}

\begin{verbatim}
## The following object is masked from 'package:stats':
## 
##     filter
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(}\StringTok{"MVTests"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Attaching package: 'MVTests'
\end{verbatim}

\begin{verbatim}
## The following object is masked from 'package:datasets':
## 
##     iris
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(}\StringTok{"heplots"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Loading required package: broom
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# salary data example}
\NormalTok{salaryData }\OtherTok{=}\NormalTok{ readxl}\SpecialCharTok{::}\FunctionTok{read\_xlsx}\NormalTok{(}\StringTok{"example.xlsx"}\NormalTok{)}
\NormalTok{DV }\OtherTok{=} \FunctionTok{cbind}\NormalTok{(salaryData}\SpecialCharTok{$}\NormalTok{salary, salaryData}\SpecialCharTok{$}\NormalTok{salbegin)}
\NormalTok{jobcat }\OtherTok{=} \FunctionTok{as.factor}\NormalTok{(salaryData}\SpecialCharTok{$}\NormalTok{jobcat)}
\NormalTok{gender }\OtherTok{=} \FunctionTok{as.factor}\NormalTok{(salaryData}\SpecialCharTok{$}\NormalTok{gender)}
\end{Highlighting}
\end{Shaded}

\hypertarget{visualization-1}{%
\subsection{visualization}\label{visualization-1}}

Note that in order to generate the boxplot according to the
``categorical'' data, ``factor'' type should be designated in the
mapping parameter.

If not ``factored'', the scatter plot works fine with grouped effect,
but the boxplot fails.

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# eyeball the data}
\FunctionTok{library}\NormalTok{(gridExtra)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Attaching package: 'gridExtra'
\end{verbatim}

\begin{verbatim}
## The following object is masked from 'package:dplyr':
## 
##     combine
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{p1 }\OtherTok{=} \FunctionTok{ggplot}\NormalTok{(salaryData, }\FunctionTok{aes}\NormalTok{(}\AttributeTok{x =} \FunctionTok{factor}\NormalTok{(jobcat), }\AttributeTok{y =}\NormalTok{ salary)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_boxplot}\NormalTok{(}\AttributeTok{outlier.shape =} \ConstantTok{NA}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{geom\_jitter}\NormalTok{(}\AttributeTok{width =} \FloatTok{0.2}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{theme}\NormalTok{(}\AttributeTok{legend.position=}\StringTok{"top"}\NormalTok{)}
\NormalTok{p2 }\OtherTok{=} \FunctionTok{ggplot}\NormalTok{(salaryData, }\FunctionTok{aes}\NormalTok{(}\AttributeTok{x =} \FunctionTok{factor}\NormalTok{(jobcat), }\AttributeTok{y =}\NormalTok{ salbegin)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_boxplot}\NormalTok{(}\AttributeTok{outlier.shape =} \ConstantTok{NA}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{geom\_jitter}\NormalTok{(}\AttributeTok{width =} \FloatTok{0.2}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{theme}\NormalTok{(}\AttributeTok{legend.position=}\StringTok{"top"}\NormalTok{)}
\FunctionTok{grid.arrange}\NormalTok{(p1, p2, }\AttributeTok{ncol=}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{r-anova-manova_files/figure-latex/unnamed-chunk-7-1.pdf}

We can also eyeball the DVs.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(car)}
\FunctionTok{scatterplot}\NormalTok{(salary }\SpecialCharTok{\textasciitilde{}}\NormalTok{ salbegin, }\AttributeTok{data=}\NormalTok{salaryData, }\AttributeTok{ellipse=}\ConstantTok{TRUE}\NormalTok{, }
  \AttributeTok{smooth=}\FunctionTok{list}\NormalTok{(}\AttributeTok{style=}\StringTok{"lines"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\includegraphics{r-anova-manova_files/figure-latex/unnamed-chunk-8-1.pdf}

\hypertarget{assumptions}{%
\subsection{assumptions}\label{assumptions}}

\hypertarget{normality}{%
\subsubsection{normality}\label{normality}}

\hypertarget{pipe-friendly-preferred}{%
\paragraph{pipe-friendly (Preferred)}\label{pipe-friendly-preferred}}

Preferred!

Note that \texttt{mshapiro\_test()} is from the package
\texttt{rstatix}, one step further from \texttt{tidyverse}.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{|\textgreater{}} \FunctionTok{select}\NormalTok{(salary,salbegin) }\SpecialCharTok{|\textgreater{}} \FunctionTok{mshapiro\_test}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## # A tibble: 1 x 2
##   statistic  p.value
##       <dbl>    <dbl>
## 1     0.707 6.13e-28
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{|\textgreater{}} \FunctionTok{group\_by}\NormalTok{(jobcat) }\SpecialCharTok{|\textgreater{}} 
  \FunctionTok{shapiro\_test}\NormalTok{(salary,salbegin)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## # A tibble: 6 x 4
##   jobcat variable statistic        p
##    <dbl> <chr>        <dbl>    <dbl>
## 1      1 salary       0.882 4.57e-16
## 2      1 salbegin     0.931 6.10e-12
## 3      2 salary       0.818 2.86e- 4
## 4      2 salbegin     0.499 1.64e- 8
## 5      3 salary       0.929 1.72e- 4
## 6      3 salbegin     0.860 2.14e- 7
\end{verbatim}

\hypertarget{classical-not-preferred}{%
\paragraph{classical (Not Preferred)}\label{classical-not-preferred}}

Not preferred. It's not efficient to do this way.

The specific function \texttt{mshapiro.test()} requires to library the
package \texttt{mvnormtest}, which surprisingly has a single function!
Moreover, it requires the data to be transposed before being processed.

The most annoying part is that the data has to be packed purposefully.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{mvnormtest}\SpecialCharTok{::}\FunctionTok{mshapiro.test}\NormalTok{(}\FunctionTok{t}\NormalTok{(DV))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.70689, p-value < 2.2e-16
\end{verbatim}

\hypertarget{homogeniety-of-variance-covariance}{%
\subsubsection{homogeniety of variance \&
covariance}\label{homogeniety-of-variance-covariance}}

Using the Box's M-test and Levene's Test.

\hypertarget{pipe-friendly-preferred-1}{%
\paragraph{pipe-friendly (Preferred)}\label{pipe-friendly-preferred-1}}

Using \texttt{box\_m()} from the package \texttt{rstatix} to test the
homogeniety of covariance matrices, in which the first parameter is the
DVs needed to be analyzed (No group variable! This can be achieved
through \texttt{select()}), and the second is the \texttt{group}
parameter specifying the group variable.

Moreover, the function \texttt{heplots::boxM()} is helpful. Instead of
taking the formula to specify the parameters, using the ``traditional''
way to write \texttt{group\ =\ ...} works pretty well combined with
pipe-friendly environment!

The drawback of the functions is that they could not search column
variables automatically, and that should be specified manually. And more
annoying one is that, when encountering with NAs, an error will be
raised by \texttt{rstatix::box\_m()}, while that's fine with
\texttt{heplots::boxM()}.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{|\textgreater{}} \FunctionTok{select}\NormalTok{(salary,salbegin) }\SpecialCharTok{|\textgreater{}} 
  \FunctionTok{box\_m}\NormalTok{(}\AttributeTok{group =}\NormalTok{ jobcat)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## # A tibble: 1 x 4
##   statistic   p.value parameter method                                          
##       <dbl>     <dbl>     <dbl> <chr>                                           
## 1      509. 1.14e-106         6 Box's M-test for Homogeneity of Covariance Matr~
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{|\textgreater{}} \FunctionTok{select}\NormalTok{(salary,salbegin) }\SpecialCharTok{|\textgreater{}} 
\NormalTok{  heplots}\SpecialCharTok{::}\FunctionTok{boxM}\NormalTok{(}\AttributeTok{group =}\NormalTok{ jobcat)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  select(salaryData, salary, salbegin)
## Chi-Sq (approx.) = 508.68, df = 6, p-value < 2.2e-16
\end{verbatim}

Do not forget that the assumption of homogeniety of variance for each DV
is needed either! However, the traditional \texttt{levene\_test()} (or
its variants) can deal with one DV at a time, which is tedious when we
have many DVs to check.

Hopefully, we have \texttt{heplots::leveneTests(data,group)} to deal
with that issue. Combined with pipe-friendly procedures, that's
fabulous!

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{|\textgreater{}}  \FunctionTok{select}\NormalTok{(salbegin,salary) }\SpecialCharTok{|\textgreater{}} 
\NormalTok{  heplots}\SpecialCharTok{::}\FunctionTok{leveneTests}\NormalTok{(}\AttributeTok{group =}\NormalTok{ jobcat)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Levene's Tests for Homogeneity of Variance (center = median)
## 
##          df1 df2 F value    Pr(>F)    
## salbegin   2 471  68.862 < 2.2e-16 ***
## salary     2 471  51.189 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\hypertarget{classical}{%
\paragraph{classical}\label{classical}}

The most annoying part is that the data has to be packed purposefully.

An alternative way to use \texttt{heplots::boxM()} has been mentioned
above. We can renovate the efficiency to carry this test by replacing
the ``formula'' specifying-pattern with the ``group'' one.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{heplots}\SpecialCharTok{::}\FunctionTok{boxM}\NormalTok{(}\FunctionTok{cbind}\NormalTok{(salary, salbegin) }\SpecialCharTok{\textasciitilde{}} \FunctionTok{factor}\NormalTok{(jobcat), }
              \AttributeTok{data =}\NormalTok{ salaryData)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y
## Chi-Sq (approx.) = 508.68, df = 6, p-value < 2.2e-16
\end{verbatim}

\hypertarget{manova-1}{%
\subsection{MANOVA}\label{manova-1}}

Use the funtion \texttt{manova()}, which is from the baseR, no special
need to library packages.

Similarly, the most annoying part is that the DVs and IVs has to be
packed previously.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fit }\OtherTok{=} \FunctionTok{manova}\NormalTok{(DV }\SpecialCharTok{\textasciitilde{}}\NormalTok{ jobcat)}
\end{Highlighting}
\end{Shaded}

As before, simply use \texttt{summary()} of the fitted model to see the
results.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary}\NormalTok{(fit)  }\CommentTok{\#for manova}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## jobcat      2 0.67298   119.43      4    942 < 2.2e-16 ***
## Residuals 471                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\hypertarget{effect-size}{%
\subsection{effect size}\label{effect-size}}

Checking the effect size is also indispensable.

Use \texttt{eta\_squared()} from the package \texttt{effectsize}, which
contains plenty of useful functions computing effect size under various
scenarios.

Note that there's also an \texttt{eta\_squared()} from the package
\texttt{rstatix()}, but that can only deal with univariate ANOVA, it
would fail when encountering MANOVA models.

Also note that the package \texttt{effectsize} is incorporated in
\texttt{bruceR}!

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{effectsize}\SpecialCharTok{::}\FunctionTok{eta\_squared}\NormalTok{(fit)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## # Effect Size for ANOVA (Type I)
## 
## Parameter | Eta2 (partial) |       95% CI
## -----------------------------------------
## jobcat    |           0.34 | [0.30, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].
\end{verbatim}

\hypertarget{univariate-anova}{%
\subsection{univariate ANOVA}\label{univariate-anova}}

Needed if MANOVA reports significant results.

\texttt{summary.aov()} will report the univariate ANOVA respectively.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{summary.aov}\NormalTok{(fit)  }\CommentTok{\#for indDVidual anovas}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  Response 1 :
##              Df     Sum Sq    Mean Sq F value    Pr(>F)    
## jobcat        2 8.9438e+10 4.4719e+10  434.48 < 2.2e-16 ***
## Residuals   471 4.8478e+10 1.0293e+08                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##              Df     Sum Sq    Mean Sq F value    Pr(>F)    
## jobcat        2 1.7926e+10 8962772266  371.11 < 2.2e-16 ***
## Residuals   471 1.1375e+10   24151508                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\hypertarget{post-hoc-tests}{%
\subsection{post-hoc tests}\label{post-hoc-tests}}

Note that the post-hoc tests here are essentially for the univariate
ANOVA.

Note that the \texttt{games\_howell\_test()} is different from
\texttt{tukey\_hsd()}.

Note that for many DVs to take the post-hoc tests simultaneously, first
\texttt{gather()} the targeted columns, and then take the
\textbf{grouped} data for the test.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{salaryData }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{gather}\NormalTok{(}\AttributeTok{key =} \StringTok{"var"}\NormalTok{, }
                      \AttributeTok{value =} \StringTok{"value"}\NormalTok{, }
\NormalTok{                      salary, salbegin) }\SpecialCharTok{|\textgreater{}} 
  \FunctionTok{group\_by}\NormalTok{(var) }\SpecialCharTok{|\textgreater{}}
  \FunctionTok{games\_howell\_test}\NormalTok{(value}\SpecialCharTok{\textasciitilde{}}\NormalTok{jobcat) }
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
# A tibble: 6 x 9
  var      .y.   group1 group2 estimate conf.low conf.high    p.adj p.adj.signif
* <chr>    <chr> <chr>  <chr>     <dbl>    <dbl>     <dbl>    <dbl> <chr>       
1 salary   value 1      2         3100.    1746.     4455. 1.22e- 6 ****        
2 salary   value 1      3        36139.   31302.    40977. 4.25e-10 ****        
3 salary   value 2      3        33039.   28196.    37881. 4.21e-10 ****        
4 salbegin value 1      2          982.     256.     1707. 6   e- 3 **          
5 salbegin value 1      3        16162.   13539.    18784. 2.74e-10 ****        
6 salbegin value 2      3        15180.   12514.    17846. 4.77e-10 ****        
\end{verbatim}

\hypertarget{ux5bb9ux6613ux8e29ux7684ux5751}{%
\section{容易踩的坑}\label{ux5bb9ux6613ux8e29ux7684ux5751}}

\begin{pitfall}[嵌套与 MANOVA 都容易栽的几件事]
\begin{enumerate}
\item \textbf{嵌套与交叉（crossed）混淆}：当「班」嵌套在「学校」里时，A 校的「3 班」和 B 校的「3 班」不是同一回事，必须写成 \texttt{Error(school/class)}。反之如果是真正的两个交叉因子（如「教学法」× 「教师」），应该用 \texttt{Y \~{} method * teacher}。判断标准：每个 level 的子因子是否在所有上级因子里出现？
\item \textbf{MANOVA 跳过协方差齐性检查就上结果}：MANOVA 的 $p$ 值对 $\Sigma_1 = \Sigma_2 = \cdots = \Sigma_g$ 的同方差假设敏感。Box's $M$ 检验给出齐性 $p$；若 $p < 0.001$ 且各组样本量差异大，Pillai 之外的统计量（特别是 Wilks）都不可信。
\item \textbf{MANOVA 显著就停手}：MANOVA 只告诉你「至少有一个 DV 上组间有差异」。要知道**是哪个 DV** 必须再跑 univariate ANOVA，最好做 Bonferroni 校正控制族错误率。
\item \textbf{事后检验用错}：组间方差齐时用 Tukey HSD；不齐时用 Games-Howell（\texttt{rstatix::games\_howell\_test}）。本笔记上面演示的就是 Games-Howell，因为各职位类别样本量差异很大。
\end{enumerate}
\end{pitfall}

\end{document}
