University of Southern California Dissertations and Theses

Robustness of rank-based and other robust estimators when comparing groups

(USC Thesis Other)

Robustness of rank-based and other robust estimators when comparing groups

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

ROBUSTNESS OF RANK-BASED AND OTHER
ROBUST ESTIMATORS WHEN COMPARING
GROUPS

by
Jinxia Ma

______________________________________________________

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Quantitative Psychology)

May 2016

Copyright 2016 Jinxia Ma

i

Dedication

To dreamers who dream big and fight with bravery.

ii

Acknowledgements

I would like to use this opportunity to express my deepest gratitude and appreciation to
all of the people that have given me assistance and support in making this work possible.
The first person I would like to thank is my wonderful advisor Professor Rand
Wilcox. I feel extremely fortunate to have had Professor Wilcox as my advisor during my
PhD studies at USC. Five years ago when I was in China and Professor Wilcox
interviewed me through Skype, I was so impressed by his wit, his knowledge and his
wisdom, that I looked forward to following him in research ever since. During my PhD
studies at USC, Professor Wilcox gave me tremendous freedom to look for my research
interest, and he encouraged me to try different directions. On the other hand, he has
always been there to give me guidance when I felt lost. He is not just my advisor, but also
my life guide. Whenever I started to doubt myself or lose my direction, Professor Wilcox
was always there to give me confidence and encouragement, to make me believe in
myself again, and find my path again. Thanks to him, my PhD life at USC was very
fulfilling. I’m very much honored to be Professor Wilcox’s student and this renowned
and respectable scholar will be my idol forever. I believe our cooperation in research will
continue and I look forward to learning more from him!
I am also very grateful to all of my committee members for their inspirational
comments, thoughts, guidance and understanding. Professor Richard John’s help and

iii

support has actually been there for me my entire five years here at USC. My first
semester as a PhD student at USC started with taking his class “Decision Theory”. Before
that, my knowledge was restricted to statistics, quantitative methods and probability
theories due to my math background. By taking Professor John’s class, I started to learn
how the quantitative psychology works, and to realize what is there in this area for me to
explore and to apply my quant skills to the questions in social science and psychological
studies. Due to his guidance and encouragement, my initial transition from math to quant
psychology was very smooth and successful. Professor John was also very supportive
when I applied for internship positions. He has also served as my committee member for
my second year project, my qualification exam and my dissertation. Every time as my
committee member, he has given me very valuable comments and directions to help
make my work be better. His hardworking attitude and scrupulous attention with research
are worth admiration from every one of us.
My sincere thanks also go to Professor David Walsh, who raised many insightful
questions, and gave me valuable advice when I presented this project at the proposal
stage. I can’t describe how happy I was when Professor Walsh agreed to serve as my
committee member and how much I appreciate his understanding, his help, and his full
support!
If there is an ideal image of a female scholar, in my mind, that would be
Professor Sandy Eckel. I got to know Professor Eckel when I took her statistics class at
HSC. Her style, her knowledge of statistics, and the way she presents those statistics
concepts were enchanting to me. The moment when I realized that I need an outside

iv

member for my committee the first person that jumped into my mind was her. I feel very
fortunate to have Professor Eckel on my committee. She has given me many valuable
advices and ideas which helped me improve my dissertation dramatically.
I also want to thank my parents. They do not understand English, but I still would
like to use this opportunity to express my deepest gratitude to them. Firstly I want to
thank them for giving me my life, twice. When I was diagnosed with having a severe
born disease and needed to have a high risk surgery 23 years ago, they didn’t give up on
me, they didn’t give in to fate, but they chose to fight. With their effort, in the end
something magical happened and I survived. Without their courage and fight back then,
there wouldn’t be me today completing this dissertation. In my heart, they are the true
dreamers and true fighters who dream impossible and fight unbeatable. I also want to
thank them for their continuous understanding and support. I stepped out from my
hometown Laiwu, the small city in the east of China, and then went to Shanghai, to
Holland and now to the US. Throughout my journeys, they have always been supportive
of every decision I made, and encouraged me to chase my dreams and to embrace new
adventures. I know how much they miss me and how much they wish I could be around
them and accompany them daily, but they are so wise and so understanding that they
never tried to tie me down. Without their support and understanding, I wouldn’t be able
to have made it this far. I’m not a very lucky person because I wasn’t born healthy. But at
the same time I am extremely lucky because I have them as my parents. Because of them,
I not only gained my second life but I was also able to experience the richness of life
fully.

v

I am also thankful to my friends in the Psychology Department, especially Jia Li,
Pan Wang, Zhisen Urgolites, Gin Chen, Tracey Cui and Miao Wei. With their support
and encouragement, I never felt I was alone.
Tremendous thanks also go to all the staff members in the Psychology
Department for their continuous extraordinary assistance during all these years – in
particular, Ms. Irene Takaragawa, Ms. Twyla Ponton and Mr. Gabriel Gonzalez.
Last but not least, I would like to thank my co-workers and mentors at the Walt
Disney Company - Charlie Cain, Bonnie Matosich, Wenny Tung Katzenstein and
Meghan Maro. When they knew about my working on my PhD thesis, all of them were
very proud of me, and showed their full support and excitement which gave me
tremendous encouragement. I feel very fortunate to have them at work, and I appreciate
their support and encouragement from the bottom of my heart! In a way I feel like I am
not fighting for my PhD just for myself, but also for them.
Looking back, no words are enough to express my gratitude to all of the people
who have supported and helped me. Moving ahead, I will carry on, and I look forward to
more adventures to come.

vi

Table of Contents

Dedication .................................................................................................................... i
Acknowledgements ...................................................................................................... ii
Table of Contents ......................................................................................................... vi
List of Tables ............................................................................................................... ix
List of Figures .............................................................................................................. xiv
Abbreviations ............................................................................................................... xv
Abstract ........................................................................................................................ xvi
Chapter 1 Introduction ................................................................................................. 1
1.1 Traditional OLS regression and its limitations .......................................... 1
1.2 Robust regression methods considered ...................................................... 3
1.3 Organization of the chapters ...................................................................... 4

Chapter 2 Description of the robust methods considered ............................................ 6
2.1 Description of the rank-based estimator .................................................... 6
2.2 Description of the Theil-Sen estimator ...................................................... 7
2.3 Description of the Modified Theil-Sen estimator ...................................... 10
2.4 Description of the MM-estimator .............................................................. 13
2.5 Description of the Least Trimmed Squares estimator................................ 15

Chapter 3 Efficiency of the robust regression methods ............................................... 16
3.1 Relative Efficiency..................................................................................... 16
3.2 Simulation Design ...................................................................................... 17
3.3 Simulation Results ..................................................................................... 19
3.4 Discussion .................................................................................................. 23

Chapter 4 Comparing regression slopes with robust methods ..................................... 24
4.1 Background ................................................................................................ 24
4.2 Comparing regression slopes applying bootstrap method ......................... 25
4.3 Simulation design....................................................................................... 26

vii

4.3.1 Simulation part 1 – test of Type I error probability .................... 27
4.3.2 Simulation part 2 – power analysis ............................................. 30
4.4 Simulation results....................................................................................... 31
4.4.1 Simulation results of the actual Type I error rates ...................... 31
4.4.1.1 Type I error rates with small sample sizes ..................... 32
4.4.1.2 Type I error rates with increased sample sizes .............. 34
4.4.2 Power analysis results ................................................................. 37
4.4.2.1 Power analysis results with small sample sizes ............. 37
4.4.2.2 Power analysis results with increased sample sizes ....... 39
4.4.2.3 Power analysis results under heteroscedasticity ............ 40
4.5 Discussion .................................................................................................. 42

Chapter 5 ANCOVA: robustness of rank-based and other regression estimators
as well as non-parametric estimator on design points ................................................. 44
5.1 Background ................................................................................................ 44
5.2 ANCOVA with the goal of comparing groups in terms of some
measure of location ........................................................................................ 46
5.3 Methods based on specific design points ................................................... 47
5.3.1 A non-parametric method ........................................................... 47
5.3.2 Comparison at design points with robust regression methods .... 49
5.3.3 Comparison at multiple points – Hochberg’s method ................ 50
5.3.4 Choose design points in a simulation study ................................ 51
5.4 Simulation design....................................................................................... 52
5.4.1 Simulation part 1 – test of Type I error probability .................... 53
5.4.2 Simulation part 2 – power analysis ............................................. 55
5.5 Simulation results....................................................................................... 56
5.5.1 Simulation results of the actual Type I error rates ...................... 56
5.5.1.1 Type I error rates with small sample sizes ..................... 56
5.5.1.2 Type I error rates with increased sample sizes .............. 58
5.5.2 Power analysis results ................................................................. 62
5.5.2.1 Power analysis with small sample sizes......................... 62
5.5.2.2 Power analysis with increased sample sizes .................. 64
5.5.2.3 Power analysis under heteroscedasticity ........................ 65
5.6 Concluding Remarks .................................................................................. 67

Chapter 6 Impact of intervention on the association between cortisol awakening
response and wellbeing measures ................................................................................ 69
6.1 Study background ...................................................................................... 69
6.2 The participants .......................................................................................... 70
6.3 The intervention ......................................................................................... 70
6.4 Measures and variables .............................................................................. 71
6.4.1 Measure of CAR ......................................................................... 71
6.4.2 Key measures of wellbeing ......................................................... 72
6.5 Analysis and results ................................................................................... 73

viii

6.5.1 Analysis of comparing the association of CAR and wellbeing
before and after the intervention .......................................................... 73
6.5.2 Results of comparing the association of CAR and wellbeing
before and after the intervention .......................................................... 74
6.6 Discussion .................................................................................................. 80

Bibliography ................................................................................................................ 83
Appendix A: R Functions ............................................................................................ 93
A.1 R function reg2ci ....................................................................................... 93
A.2 R function ancts ........................................................................................ 95
A.3 R function ancova ..................................................................................... 96
Appendix B: 20%-Trimmed Mean .............................................................................. 97
Appendix C: True power for comparing groups at design points with parallel
slopes............................................................................................................................ 98
Appendix D: Simulation results when X and ε have different distributions ............... 102
D.1 Relative efficiency .................................................................................... 102
D.2 Compare slopes between two groups ........................................................ 105
D.3 Compare two groups at design points ....................................................... 108

ix

List of Tables

Table 3-1: Relative efficiency with sample size n=20 and
under homoscedasticity……….……………...…………………………..20

Table 3-2: Relative efficiency with sample size n=20 and
under heteroscedasticity……….……………...………………………….20

Table 3-3: Relative efficiency with sample size n=30 and
under homoscedasticity……….……………...…………………………..21

Table 3-4: Relative efficiency with sample size n=30 and
under heteroscedasticity……….……………...………………………….22

Table 3-5: Relative efficiency with sample size n=40 and
under homoscedasticity……….……………...…………………………..22

Table 3-6: Relative efficiency with sample size n=40 and
under heteroscedasticity……….……………...………………………….23

Table 4-1: Some properties of the g-and-h distribution……………………………...29

Table 4-2: Summary of the simulation design………..……………………………...30

Table 4-3: Actual Type I error rates under homoscedasticity
(𝑛 1
=20, 𝑛 2
=20)………………………………….…………………...32

Table 4-4: Actual Type I error rates under heteroscedasticity
(𝑛 1
=20, 𝑛 2
=20)………………………………….…………………...33

Table 4-5: Actual Type I error rates under homoscedasticity
(𝑛 1
=20, 𝑛 2
=40)………………………………….…………………...35

Table 4-6: Actual Type I error rates under heteroscedasticity

x

(𝑛 1
=20, 𝑛 2
=40)………………………………….…………………...35

Table 4-7: Actual Type I error rates under homoscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...36

Table 4-8: Actual Type I error rates under heteroscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...37

Table 4-9: Power analysis results under homoscedasticity
(𝑛 1
=20, 𝑛 2
=20)………………………………….…………………...38

Table 4-10: Power analysis results under heteroscedasticity
(𝑛 1
=20, 𝑛 2
=20)………………………………….…………………...38

Table 4-11: Power analysis results under homoscedasticity
(𝑛 1
=20, 𝑛 2
=40)………………………………….…………………...39

Table 4-12: Power analysis results under heteroscedasticity
(𝑛 1
=20, 𝑛 2
=40)………………………………….…………………...40

Table 4-13: Power analysis results under homoscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...41

Table 4-14: Power analysis results under heteroscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...41

Table 5-1: Summary of the simulation design………..……………………………...55

Table 5-2: Actual Type I error rates under homoscedasticity
(𝑛 1
=30, 𝑛 2
=30)………………………………….…………………...57

Table 5-3: Actual Type I error rates under heteroscedasticity
(𝑛 1
=30, 𝑛 2
=30)………………………………….…………………...58

Table 5-4: Actual Type I error rates under homoscedasticity
(𝑛 1
=30, 𝑛 2
=40)………………………………….…………………...59

Table 5-5: Actual Type I error rates under heteroscedasticity
(𝑛 1
=30, 𝑛 2
=40)………………………………….…………………...60

Table 5-6: Actual Type I error rates under homoscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...61

Table 5-7: Actual Type I error rates under heteroscedasticity

xi

(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...61

Table 5-8: Power analysis results under homoscedasticity
(𝑛 1
=30, 𝑛 2
=30)………………………………….…………………...63

Table 5-9: Power analysis results under heteroscedasticity
(𝑛 1
=30, 𝑛 2
=30)………………………………….…………………...63

Table 5-10: Power analysis results under homoscedasticity
(𝑛 1
=30, 𝑛 2
=40)………………………………….…………………...64

Table 5-11: Power analysis results under heteroscedasticity
(𝑛 1
=30, 𝑛 2
=40)………………………………….…………………...65

Table 5-12: Power analysis results under homoscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...66

Table 5-13: Power analysis results under heteroscedasticity
(𝑛 1
=40, 𝑛 2
=40)………………………………….…………………...66

Table 6-1: Compare association of CAR and MAPA
before and after intervention (outliers in IV
included)……………………………………………...……………..…...77

Table 6-2: Compare association of CAR and MAPA
before and after intervention (outliers in IV
eliminated)…………………………………………...……………..……78

Table 6-3: Compare association of CAR and LSIZ
before and after intervention (outliers in IV
included)……………………………………………...……………..…...78

Table 6-4: Compare association of CAR and LSIZ
before and after intervention (outliers in IV
eliminated)……………………………………………...…………..……79

Table 6-5: Compare association of CAR and PC
before and after intervention (outliers in IV
included)……………………………………………...……………..…...79

Table 6-6: Compare association of CAR and PC
before and after intervention (outliers in IV
eliminated)……………………………………………...…………..……80

xii

Table C-1: True Power (with parallel slopes;
𝑛 1
=30,𝑛 2
=30; with homoscedasticity)……………...…………..…...99

Table C-2: True Power (with parallel slopes;
𝑛 1
=30,𝑛 2
=30; with heteroscedasticity)……………...…………..…..99

Table C-3: True Power (with parallel slopes;
𝑛 1
=30,𝑛 2
=40; with homoscedasticity)……………...…………..…...100

Table C-4: True Power (with parallel slopes;
𝑛 1
=30,𝑛 2
=40; with heteroscedasticity)……………...…………..…..100

Table C-5: True Power (with parallel slopes;
𝑛 1
=40,𝑛 2
=40; with homoscedasticity)……………...……………….101

Table C-6: True Power (with parallel slopes;
𝑛 1
=40,𝑛 2
=40; with heteroscedasticity)……………...…………..…..101

Table D-1: Relative Efficiency – RE = OLS / Robust Methods
( n =20; with homoscedasticity)……………...……………………...…..102

Table D-2: Relative Efficiency – RE = OLS / Robust Methods
( n =20; with heteroscedasticity)…………......……………………...…..103

Table D-3: Relative Efficiency – RE = OLS / Robust Methods
( n =30; with homoscedasticity)……………...……………………...…..103

Table D-4: Relative Efficiency – RE = OLS / Robust Methods
( n =30; with heteroscedasticity)…………......……………………...…..103

Table D-5: Relative Efficiency – RE = OLS / Robust Methods
( n =40; with homoscedasticity)……………...……………………...…..104

Table D-6: Relative Efficiency – RE = OLS / Robust Methods
( n =40; with heteroscedasticity)…………......……………………...…..104

Table D-7: Actual Type I Error Rates
(𝑛 1
=20,𝑛 2
=20; with homoscedasticity)…………......…………...…...105

Table D-8: Actual Type I Error Rates
(𝑛 1
=20,𝑛 2
=20; with heteroscedasticity)…………......…………...…..105

xiii

Table D-9: Actual Type I Error Rates
(𝑛 1
=20,𝑛 2
=40; with homoscedasticity)…………......…………...…..106

Table D-10: Actual Type I Error Rates
(𝑛 1
=20,𝑛 2
=40; with heteroscedasticity)…………......…………...….106

Table D-11: Actual Type I Error Rates
(𝑛 1
=40,𝑛 2
=40; with homoscedasticity)…………......…………...…..106

Table D-12: Actual Type I Error Rates
(𝑛 1
=40,𝑛 2
=40; with heteroscedasticity)…………......…………...….107

Table D-13: Actual Type I Error Rates
(𝑛 1
=30,𝑛 2
=30; with homoscedasticity)…………......…………...…..108

Table D-14: Actual Type I Error Rates
(𝑛 1
=30,𝑛 2
=30; with heteroscedasticity)…………......…………...….108

Table D-15: Actual Type I Error Rates
(𝑛 1
=30,𝑛 2
=40; with homoscedasticity)…………......…………...…..109

Table D-16: Actual Type I Error Rates
(𝑛 1
=30,𝑛 2
=40; with heteroscedasticity)…………......…………...….109

Table D-17: Actual Type I Error Rates
(𝑛 1
=40,𝑛 2
=40; with homoscedasticity)…………......…………...…..109

Table D-18: Actual Type I Error Rates
(𝑛 1
=40,𝑛 2
=40; with heteroscedasticity)…………......…………...….110

xiv

List of Figures

Figure 3-1: g-and-h distribution……………………………………………………...18

Figure 6-1: Association between CAR and MAPA
before and after intervention estimated
with the Rank-Based method (with
outliers in IV eliminated…………………………………………………75

Figure 6-2: Association between CAR and LSIZ
before and after intervention estimated
with the Rank-Based method (with
outliers in IV eliminated…………………………………………………76

Figure 6-3: Association between CAR and PC
before and after intervention estimated
with the Rank-Based method (with
outliers in IV eliminated…………………………………………………77

xv

Abbreviations

OLS: Ordinary Least Squares
LTS: Least Trimmed Squares
RE: Relative Efficiency
VP: Variance Pattern
IV: Independent Variable
DV: Dependent Variable
HPA: Hypothalamus-Pituitary-Adrenal
CAR: Cortisol Awakening Response
MAPA: Meaningful Activity Participation Assessment
PC: Perceived Control
LSIZ: Life Satisfaction Index-Z
HOM: Homoscedasticity
HET: Heteroscedasticity

xvi

Abstract

This thesis studied robust regression methods, especially rank-based method and
its performance compared to other robust regression methods (the Theil-Sen method, the
modified Theil-Sen method, the Least Trimmed Squares method and the MM-Estimator)
and the traditional OLS method. Three simulation studies were conducted to determine
the efficiency of the methods when estimating the regression slopes, the robustness of the
methods when comparing slopes between two groups and the robustness of the methods
when comparing slopes at design points between two groups.
Generally speaking, in terms of efficiency, the rank-based method and the MM-
Estimator have better efficiency than the other robust regression methods under
homoscedasticity. But the Theil-Sen and the modified Theil-Sen method have better
efficiency under heteroscedasticity.
As far as comparing regression lines, the results indicate that three robust methods
(the Rank-Based method, the Theil-Sen method and the modified Theil-Sen method)
perform well in the simulation studies, even when there is heteroscedasticity. Another
robust method (the MM-Estimator) performs satisfactorily as well when sample size is
relatively large. These four methods studied here are recommended for general use.

xvii

When it comes to comparing the slopes at design points, no robust regression
method performs well in Type I error probability control when the sample size is
relatively small. When the sample size is relatively large, the Theil-Sen and the modified
Theil-Sen method perform well but tend to be conservative. Non-Parametric method
provides good control over Type I error probability across the distributions and variance
patterns simulated; however, its power is low across simulated distributions.
Lastly, the robust regression methods being considered and the OLS method were
applied to a real data set to compare the association between CAR and wellbeing
measures before and after the intervention. Results were reported and a note was made
that choosing the right method is the way to help us avoid misleading results and drawing
wrong conclusions.

1

Chapter 1

Introduction

1.1 Traditional OLS Method and Its Limitations

To estimate the association between two variables, the traditional method is to run an
OLS (ordinary least squares) regression. However, the OLS estimator can give us very
misleading results if the assumptions required are violated. In more detail, suppose (yi,
xi1, · · · xip), i = 1, · · · , n, are n vectors of observations randomly sampled from some
(p+1)-variate distribution, where (xi1, · · · xip) is a vector of predictor values. In
regression, the most common assumption is that
yi = β0 + β1xi1 +· · ·+ βpxip+ εi, (Eq.1-1)
where ( β0, · · · , βp) are unknown parameters, i = 1, · · · , n, the εi are independent random
variables with normal distribution, and 𝜀 𝑖 is independent of 𝑥 𝑖𝑘
for each k, k = 1, · · · , p.
Let 𝑏 𝑗 be any estimate of βj, j = 0, 1, · · · , p, and let
𝑦̂
𝑖 =𝑏 0
+𝑏 1
𝑥 𝑖 1
+⋯+𝑏 𝑝 𝑥 𝑖𝑝
.

2

This model implies that the conditional mean of 𝑦 𝑖 , given (𝑥 𝑖 1
,…,𝑥 𝑖𝑝
), is 𝛽 0
+
∑𝛽 𝑘 𝑥 𝑖𝑘
, a linear combination of the predictors. Equation (1-1) is a homoscedastic model,
meaning that the 𝜀 𝑖 have a common variance. If the error term, 𝜀 𝑖 , has variance 𝜎 𝑖 2
, and
𝜎 𝑖 2
≠𝜎 𝑗 2
, for some 𝑖 ≠𝑗 , the model is said to be heteroscedastic.
The OLS estimates are the b values that minimize
∑(𝑦 𝑖 −𝑏 1
𝑥 𝑖 1
−⋯𝑏 𝑝 𝑥 𝑖𝑝
−𝑏 0
)
2
, (Eq.1-2)
the sum of squared residuals.
In order to test hypotheses or compute confidence intervals, typically the
homoscedastic model given by Equation (1-1) is assumed with the additional assumption
that ε has a normal distribution with mean zero. Even when ε is normal, but
heteroscedastic, problems with computing confidence intervals arise. More specifically,
when one or more of these assumptions are violated, the OLS estimation may yield poor
probability coverage.
Another problem with OLS is that it can have large standard errors which leads to
high inefficiency, and this can result in relatively low power. For example, the variance
of the error term becomes inflated if sampling is from a heavy-tailed distribution for the
dependent variable. That is, slight departures from normality result in large increases in
the standard error of the estimated parameters. It is outliers among the dependent variable
that is the issue. Outliers among the independent variable can reduce the standard errors.
(This problem is well known and discussed by Hampel, 1973; He et al., 1990; Schrader &
Hettmansperger, 1980, among others.)

3

Yet another practical problem is that OLS has a breakdown point of only
1
𝑛 . That
is, a single point, properly placed, can cause the OLS estimator to have virtually any
value. Not only do unusual y values cause problems, outlying x values, called leverage
points, can have an inordinate influence on the estimated slopes and intercept.

1.2 Robust Regression Methods Considered in This Thesis

Robust methods are designed to deal with the practical concerns described in Section 1.1.
This justifies the needs of robust regression methods. In this thesis, such robust estimator
as rank-based estimate, Theil-Sen estimate, Modified Theil-Sen estimate, MM-estimator
and Least Trimmed Squares estimator are both studied via the simulation study and
applied to the real data analysis to compare the regression slopes across two groups. We
also compare their standard errors to get a view of the efficiency they deliver. As far as
ANCOVA, we compare the robust methods’ performances in comparing two groups at
design points via another simulation study. Detailed descriptions of these robust
estimators are provided in Chapter 2.
These five robust methods are covered in this thesis because they have been found
to perform relatively well in past studies but it is unknown how they compare to the rank-
based method. Accordingly, the goal in this dissertation is to fill this gap.
For more information about robust regression, see Belsley, Kuh, and Welsch
(1980), Birkes and Dodge (1993), Carroll and Ruppert (1988), Cook and Weisberg

4

(1992), Fox (1999), Hampel et al. (1986), Hettmansperger (1984), Hettmansperger and
McKean (1998), Huber (1981), Li (1985), Long and Ervin (2000), Mammen (1993),
Maronna, Martin, and Yohai (2006), Montgomery and Peck (1992), Nanayakkara and
Cressie (1991), Rousseeuw and Leroy (1987), Staudte and Sheather (1990) and Wilcox
(1996c).

1.3 Organization of the Chapters

The dissertation starts with introducing the traditional OLS regression method and
addressing its limitations. There naturally calls for the great need in robust regression
methods. Then five robust regression methods covered in this thesis are introduced and
described in details in Chapter 2.
Chapter 3 reports a simulation test on the five robust methods to provide us a
view of their efficiency when estimating regression slopes. Simulation test continues in
Chapter 4 and Chapter 5. In Chapter 4, robustness of the 6 methods (5 robust methods
and the OLS method) are tested and compared when comparing regression slopes
between two groups. The simulation test of Chapter 5 is to compare the regression slopes
at design points. The methods compared include the 5 robust methods introduced in
Chapter 2, the OLS method and a non-parametric method.
Lastly, in Chapter 6, the robust regression methods are applied to a real data set.
The goal of the study is to figure out the impact of an intervention on the association of

5

cortisol awakening response and wellbeing measures through comparing the association
before and after the intervention. The robust regression methods as well as the traditional
OLS method are applied to conduct this analysis and results are reported.

6

Chapter 2

Description of the Robust Methods That Were Considered

2.1 Description of the Rank-Based Estimator

The rank-based estimator we are testing is from Jaeckel's study in 1972. Traditional OLS
method is used to minimize equation (1-2), the sum of squared residuals. In other words,
OLS method is used to minimize the Euclidian distance between the observed and
estimated values. To obtain the rank-based estimator, Jaeckel provided a different
measure of distance. The rank based estimator of 𝛽 is defined as

𝛽 ̂
𝜑 =𝐴𝑟𝑔𝑚𝑖𝑛 ||𝑦 −𝑋𝛽 ||
𝜑 ,

Where the norm 𝜑 is defined by

7

||𝑢 ||
𝜑 =∑ 𝑎 [𝑅 (𝑦 𝑖 −𝑋 𝑖 ′
𝛽 )](𝑦 𝑖 −𝑋 𝑖 ′
𝛽 )
𝑛 𝑖 =1
,

and 𝑅 (𝑦 𝑖 −𝑋 𝑖 ′
𝛽 ) denotes the rank of 𝑦 𝑖 −𝑋 𝑖 ′
𝛽 among 𝑦 1
−𝑋 1
′
𝛽 ,…,𝑦 𝑛 −𝑋 𝑛 ′
𝛽 . The rank
scores are given by 𝑎 (𝑖 )=𝜑 (
𝑖 𝑛 +1
) where 𝜑 (𝑢 ) is a non-decreasing and squared-
integrable function defined on the interval 0<𝑢 <1.
Without loss of generality we assume that the scores are standardized. The score
function that was used in this study is the Wilcoxon linear score function φ(u)=
√12(u−
1
2
) . The only assumption of the error distributions here is the finite Fisher
information condition. Under this assumption, β
̂
φ
is distributed approximately
N(β,τ
φ
2
(X
′
X)
−1
) , where τ
φ
is a scale parameter (Jaeckel, 1972; Hettmansperger &
McKean, 2011).
In this study, we used the function Rfit which implemented the estimation of the
parameter τ
φ
developed by Koul et al. (1987).

2.2 Description of the Theil-Sen Estimator

The Theil-Sen estimator was first proposed by Theil (1959) and later extended by Sen
(1968) that is restricted to the case of a single predictor (𝑝 =1) .Then various extensions
to 𝑝 >1 are discussed.

8

To get the estimated regression slope of with the Theil-Sen estimator, the basic
idea is to draw a line through any two distinctive data points and get their regression
slope for these two points. Then use the median of all of the obtained regression slopes as
the estimated slope for this data set. The estimated intercept is the difference between the
median of the dependent variable y and the multiplication of the estimated slope and the
median of the independent variable X.
To elaborate, temporarily focusing on 𝑝 =1, one view of the Theil-Sen estimator
is that it attempts to find a value for the slope that makes Kendall's correlation 𝜏 , between
𝑦 𝑖 −𝑏 𝑥 𝑖 and 𝑥 𝑖 , (approximately) equal to zero. This can be seen to be tantamount to the
following method. For any 𝑖 ,𝑖 ′
, for which 𝑥 𝑖 ≠𝑥 𝑖 ′, let

𝑆 𝑖 𝑖 ′ =
𝑌 𝑖 −𝑌 𝑖 ′
𝑋 𝑖 −𝑋 𝑖 ′
.

The Theil-Sen estimate of the slope is 𝑏 1𝑡𝑠
, the median of all the slopes
represented by 𝑆 𝑖 𝑖 ′. The intercept is estimated with

𝑀 𝑦 −𝑏 1𝑡𝑠
𝑀 𝑥 ,

where 𝑀 𝑦 and 𝑀 𝑥 are the usual sample medians of the y and x values, respectively.
Sen (1968) derived the asymptotic standard error of the slope estimator. Dietz
(1978) showed that the Theil-Sen estimator has an asymptotic breakdown point of 0.293.

9

For results on its small-sample efficiency, see Dietz (1989), Talwar (1991) and Wilcox
(1998a, 1998b). Sievers (1978) and Scholz (1978) proposed a generalization of the Theil-
Sen estimator by attaching weights to the pairwise slopes, 𝑆 𝑖 𝑖 ′, but in terms of efficiency
it seems to offer little or no advantage, and in some cases its bias is considerably larger,
so it is not described here. For results when the median of the slopes is replaced by a
trimmed mean, see Luh and Guo (2000).
It is noted that the Theil-Sen estimator has close similarities to a regression
estimator studied by Maronna and Yohai (1993, Section 3.3). This alternative estimator
belongs to what are called projection estimators; see "Introduction to Robust Estimation
and Hypothesis Testing" by Wilcox R. Martin, Yohai, and Zamar (1989) proved that this
alternative estimate of the slope is minimal in the class of all regression equivariant
estimates if the error term has a symmetric distribution.
There are at least three general ways the Theil-Sen estimator might be extended to
two or more predictors (cf. Hussain & Sprent, 1983). One of them is called method TS, is
to apply the back-fitting, Gauss-Seidel method and its properties, see e.g., Dahlquist &
Bjorck, 1974; Golub & van Loan, 1983.) For the problem at hand, the method is applied
as follows:
1. Set 𝑘 =0 and choose an initial estimate for 𝛽 𝑗 , say 𝑏 𝑗 (0)
, 𝑗 =0,…,𝑝 . Here, the
initial estimate is taken to be the Theil-Sen estimate of the slope based on the
jth regressor only. That is, simply ignore the other predictors to obtain an
initial estimate of 𝛽 𝑗 . The initial estimate of 𝛽 0
is the median of

10

𝑦 𝑖 −∑ 𝑏 𝑗 (0)
𝑥 𝑖𝑗
𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 .

2. Increment k and take the kth estimate of 𝛽 𝑗 (𝑗 =1,…,𝑝 ) to be 𝑏 𝑗 (𝑘 )
, the Theil-
Sen estimate of the slope based on the regression estimate of 𝑥 𝑖𝑗
with

𝑟 𝑖 =𝑦 𝑖 −𝑏 0
(𝑘 −1)
−∑ 𝑏 𝑗 (𝑘 −1)
𝑥 𝑖𝑙
𝑝 𝑙 =1,𝑙 ≠𝑗 .

The updated estimate of the intercept, 𝑏 0
(𝑘 )
, is the median of

𝑦 𝑖 −∑ 𝑏 𝑗 (𝑘 )
𝑥 𝑖𝑗
𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 .

3. Repeat step 2 until convergence.

2.3 Description of the Modified Theil-Sen Estimate

When the dependent variable is continuous, the Theil-Sen estimator enjoys good
theoretical properties and it performs well in simulations in terms of power and Type I
error probabilities when testing hypotheses about the slope (e.g., Wilcox, 2012b). Its
mean squared error and small-sample efficiency compare well to the least squares
estimator as well as other robust estimators that have been derived (Dietz, 1987; Wilcox,

11

1998). Dietz (1989) established that its asymptotic breakdown point is approximately .29.
Other asymptotic properties have been studied by Wang (2005) and Peng et al. (2008).
Akritas et al. (1995) applied it to astronomical data and Fernandes and Leblanc (2005) to
remote sensing.
Although the bulk of the results on the Theil-Sen estimator deal with situations
where the dependent variable is continuous, an exception is the paper by Peng et al.
(2008) that includes results when dealing with a discontinuous error term. They show that
when the distribution of the error term is discontinuous, the Theil-Sen estimator can be
super inefficient. They establish that even in the continuous case, the slope estimator may
or may not be asymptotically normal. Peng et al. also establish the strong consistency and
the asymptotic distribution of the Theil-Sen estimator for a general error distribution.
Currently, a basic percentile bootstrap seems best when testing hypotheses about the
slope and intercept, which has been found to perform well even when the error term is
heteroscedastic (e.g., Wilcox, 2012b).
The Theil--Sen estimate of the slope is the usual sample median based on all of
the slopes associated with any two distinct points. There is, however, a practical concern:
when there are tied values, the sample median is not necessarily asymptotically normal.
Rather, as sample size increases, the cardinality of its sample can decrease, which in turn
creates concerns about the more obvious methods for testing hypotheses.
Recent results on comparing quantiles (Wilcox et al., 2013) suggest a
modification that might deal with the concerns previously indicated: replace the usual

12

sample median with the Harrell and Davis (1982) estimate of the median, which uses a
weighted average of all the order statistics.
To describe the computational details, let (𝑌 1
,𝑋 1
),…,(𝑌 𝑛 ,𝑋 𝑛 ) be a random
sample from some unknown bivariate distribution. Assuming that 𝑋 𝑗 ≠𝑋 𝑘 for any 𝑗 <𝑘 ,
let

𝑏 𝑗𝑘
=
𝑌 𝑗 −𝑌 𝑘 𝑋 𝑗 −𝑋 𝑘 , 1≤𝑗 <𝑘 ≤𝑛 .

The Theil-Sen estimate of the slope, 𝛽̂
1
, is taken to be the usual sample median
based on the 𝑏 𝑗𝑘
values. The intercept is typically estimated with 𝛽 ̂
0
=𝑀 𝑦 −𝛽 ̂
1
𝑀 𝑥 ,
where 𝑀 𝑦 is the usual sample median based on 𝑌 1
,…,𝑌 𝑛 . This will be called the TS
estimator henceforth.
For notational convenience, let 𝑍 1
,…,𝑍 𝑙 denote the 𝑏 𝑗𝑘
values, where 𝑙 =(𝑛 2
−
𝑛 )/2. Let U be a random variable having a beta distribution with parameters 𝑎 =
(𝑙 +1)𝑞 and 𝑏 =(𝑙 +1)(1−𝑞 ) , 0<𝑞 <1. Let

𝑊 𝐼 =𝑃 (
𝑖 −1
𝑙 ≤𝑈 ≤
𝑖 𝑙 ) .

Let 𝑍 (1)
≤⋯≤𝑍 (𝑙 )
denote the 𝑍 1
,…,𝑍 𝑙 values written in ascending order. The
Harrell and Davis (1982) estimate of the qth quantile is

13

𝜃̂
𝑞 =∑𝑊 𝑖 𝑍 (𝑖 )
.

Consequently, estimate the slope with 𝛽̃
1
=𝜃̃
.5
. The intercept is estimated with
the Harrell-Davis estimate of the median based on 𝑌 1
−𝛽̃
1
𝑋 1
,…,𝑌 𝑛 −𝛽̃
𝑛 𝑋 𝑛 .
This will be called the HD estimator or the Modified Theil-Sen method.
So the strategy is to avoid the problem associated with the usual sample median
by using a quantile estimator that results in a sampling distribution that in general does
not have tied values. Because the Harrell-Davis estimator uses all of the order statistics,
the expectation is that in general it accomplishes this goal.

2.4 Description of the MM-Estimator

This method is to compute MM regression estimate derived by Yohai (1987). The MM-
estimator derived by Yohai (1987) has certain similarities to the generalized M-
estimators. It has the highest possible breakdown point, 0.5, and high efficiency under
normality. The parameters are estimated by solving an equation similar to

∑ (𝜔 𝑖 𝛹 (
𝑟 𝑖 𝜔 𝑖 𝜏̂
)𝑥 𝑖𝑗
)
𝑛 𝑖 =1
=0 ,

14

with 𝛹 taken to be some redescending function. A popular choice is Tukey's biweight,
which will be used here unless stated otherwise. That is, the regression parameters are
estimated by determining the solution to 𝑝 +1 equations

∑ 𝛹 (
𝑟 𝑖 𝜏̂
)𝑥 𝑖𝑗
𝑛 𝑖 =1
=0 .

𝑗 =0,…,𝑝 , where again 𝜏 ̂ is a robust measure of variation based on the residuals and

𝛹 (𝑟 𝑖 ;𝑐 )=
𝑟 𝑖 𝜏̂
[(
𝑟 𝑖 𝑐 𝜏̂
)
2
−1]
2
, if |
𝑟 𝑖 𝜏̂
|≤𝑐 ;

otherwise 𝛹 (𝑟 𝑖 ;𝑐 )=0.
The choice c=4.685 leads to an MM-estimator with 95% efficiency compared
to the least squares estimator and is the default value used here. The ratio p/n is relevant
to the efficiency of this estimator; see Maronna and Yohai (2010) for details. In addition
to having excellent theoretical properties, the small-sample efficiency of the MM-
estimator appears to compare well with other robust estimators, but like several other
robust estimators it can be sensitive to what is called contamination bias. Also, situations
are encountered where the iterative estimation scheme used to compute the MM-
estimator does not converge.

15

2.5 Description of the Least Trimmed Squares Estimator

Rousseeuw's (1984) least trimmed squares (LTS) estimator is based on minimizing

∑ 𝑟 (𝑖 )
2 ℎ
𝑖 =1
,

where 𝑟 (1)
2
≤⋯≤𝑟 (ℎ)
2
are the squared residuals written in ascending order. Instead of
sum up all of the squared residuals as the OLS method does, the LTS estimator only sums
up some of the squared residuals but trims off the rest. The key is to decide which ones
are to be kept and which ones are to be trimmed off. The idea is to sort all of the squared
residuals in ascending order, and only keep the first
𝑛 2
+
𝑝 +1
2
squared residuals. (With
ℎ=
𝑛 2
+1, the same breakdown point as LMS (Least Median of Squares), which is
approximately 0.5 is achieved.) However, ℎ=
𝑛 2
+
𝑝 +1
2
is often used to maintain
regression equivariance. LTS has a relatively low asymptotic efficiency (Croux,
Rousseeuw, & Hossjer, 1994), but it seems to have practical value. For example, it plays
a role in the asymptotically efficient M-estimator.

16

Chapter 3

Simulation Test of Efficiency of the Robust Regression Methods

3.1 Relative Efficiency

The goal of this section is to test the efficiency of the robust methods described in
Chapter 2 through computing their relative efficiency via a simulation test. The relative
efficiency (RE) is defined as follows: Suppose 𝑇 (1)
and 𝑇 (2)
are two estimators for which
(𝑇 (1)
)=𝐸 (𝑇 (2)
)=𝜇 , then relative efficiency (RE) is

𝑅𝐸 =
𝑉𝐴𝑅 (𝑇 (1)
)
𝑉𝐴𝑅 (𝑇 (2)
)

and it is called the relative efficiency of estimator 2 to estimator 1. When the estimators
are not unbiased, it is standard to compute

17

𝑅𝐸 =
𝑀𝑆𝐸 (𝑇 (1)
)
𝑀𝑆𝐸 (𝑇 (2)
)

for the relative efficiency. In either case 𝑅𝐸 <1 means estimator 1 is preferred. Or in
other words, estimator 2 is inefficient relative to estimator 1 in this sense.

3.2 Simulation Design

Data are generated from the model

𝑌 =𝑋 +𝐻 (𝑋 )𝜀 ,
where H is a function of X used to model heteroscedasticiy. The independent variable X
and the error term 𝜖 are generated from the g-and-h distributions (Hoaglin, 1985), which
contain the standard normal distribution as a special case. If Z has a standard normal
distribution, then

𝑊 =
{

exp(𝑔𝑍 )−1
𝑔 exp(
ℎ𝑍 2
2
), 𝑖𝑓 𝑔 >0
𝑍𝑒𝑥𝑝 (
ℎ𝑍 2
2
), 𝑖𝑓 𝑔 =0

has a g-and-h distribution.

18

When 𝑔 =0 and ℎ=0, W has a standard normal distribution. Skewness and
heavy-tailedness of the g-and-h distribution are determined by the values of g and h,
respectively. As the value of g increases, the distribution becomes more skewed. As the
value of h increases, the distribution becomes more heavy-tailed.
Four types of distributions are considered for X: standard normal (g=0, h=0),
asymmetric light-tailed (g=0.2, h=0), symmetric heavy-tailed (g=0, h=0.2), and
asymmetric heavy-tailed (g=0.2, h=0.2). More properties of the g-and-h distribution can
be found in Table 4-1. The error term 𝜖 is also randomly generated based on one of these
four g-and-h distributions. Additional properties of the g-and-h distribution are
summarized by Hoaglin (1985). To simulate homosecedasticity, we let 𝐻 (𝑋 )=1. To
simulate heteroscedasticity, we let 𝐻 (𝑋 )=|𝑋 |+1.

Figure 3-1. g-and-h distributions

-4 -2 0 2 4
0.0 0.1 0.2 0.3 0.4
Standard Normal
-4 -2 0 2 4 6
0.0 0.1 0.2 0.3 0.4
Asymmetric Light-Tailed,g=0.2,h=0
-4 -2 0 2 4
0.0 0.1 0.2 0.3 0.4
Symmetric Heavy-Tailed,g=0,h=0.2
-4 -2 0 2 4 6 8
0.0 0.1 0.2 0.3 0.4 Asymmetric Heavy-Tailed,g=0.2,h=0.2

19

The sample sizes considered are: 𝑛 =20,𝑛 =30,and 𝑛 =40 respectivly. The
relative efficiency is computed as:

𝑅𝐸 =
MSE(OLS Estimator)
MSE(Robust Estimator)

Simulation is based on 5,000 replications.

3.3 Simulation Results

Table 3-1 to Table 3-6 report the relative efficiency results of the robust methods
described in Chapter 2 to the OLS estimator in estimating regression slopes. When
sample size is small, see Table 3-1 with sample size n=20 under homoscedasticity, the
Rank Based estimator and the MM-estimator have relatively better efficiency, followed
by the Theil-Sen and the Modified Theil-Sen estimator, and then followed by the LTS
estimator. There is a big gap between the efficiency of the LTS estimator and the other
four robust estimators.
When the sample size is still small, n=20, but the variance pattern is
heteroscedasticity (see Table 3-2), the Theil-Sen and the modified Theil-Sen estimators
have better efficiency than the Rank-Based estimator and the MM-estimator, especially

20

when the distribution is heavy-tailed. The LTS estimator again does not show
satisfactory efficiency compared to the other few robust estimators.

Table 3-1. Relative Efficiency - OLS/Robust Methods (n=𝟐𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
0.96
1.00
1.18
1.24
0.92
0.95
1.10
1.15
0.92
0.94
1.10
1.15
0.96
0.98
1.15
1.21
0.70
0.74
0.88
0.91

Table 3-2. Relative Efficiency - OLS/Robust Methods (n=𝟐𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
0.98
1.03
1.27
1.35
1.02
1.10
1.52
1.71
1.02
1.09
1.52
1.70
0.91
0.97
1.30
1.45
0.69
0.74
1.09
1.22

21

Table 3-3 and 3-4 show the efficiency results when sample size is increased to 30.
Similar to the results when sample size is smaller (n=20), when there is no
heteroscedasticity in the simulated data, Rank-Based and MM-estimator are more
preferred than Theil-Sen and modified Theil-Sen estimators. While it is the opposite
when there is heteroscedasticity, in which the opposite is true – the Theil-Sen and the
modified Theil-Sen estimators are more preferred. We see the same pattern when the
sample size is even larger (see Table 3-5 and 3-6 with n=40).

Table 3-3. Relative Efficiency - OLS/Robust Methods (n=𝟑𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
0.97
1.01
1.20
1.26
0.94
0.97
1.11
1.16
0.94
0.96
1.11
1.16
0.96
1.00
1.19
1.26
0.75
0.77
0.93
0.97

22

Table 3-4. Relative Efficiency - OLS/Robust Methods (n=𝟑𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
1.00
1.04
1.35
1.41
1.05
1.12
1.63
1.76
1.04
1.11
1.62
1.76
0.93
0.99
1.37
1.48
0.68
0.74
1.13
1.25

Table 3-5. Relative Efficiency - OLS/Robust Methods (n=𝟒𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
0.97
1.00
1.21
1.27
0.94
0.96
1.12
1.17
0.94
0.96
1.12
1.17
0.97
0.99
1.21
1.26
0.76
0.79
0.97
1.01

23

Table 3-6. Relative Efficiency - OLS/Robust Methods (n=𝟒𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
0
.2
.0
.0
.2
.2
1.01
1.06
1.44
1.55
1.06
1.13
1.73
1.93
1.06
1.13
1.73
1.93
0.93
0.99
1.44
1.58
0.67
0.72
1.19
1.33

3.4 Discussion

Generally speaking, when there is no heteroscedasticity, the Rank-Based method and the
MM-Estimator have smaller standard errors, i.e., better efficiency than the Theil-Sen and
the Modified Theil-Sen methods. When there is heteroscedasticity, it is the opposite - the
Theil-Sen and the Modified Theil-Sen methods have better efficiency than the Rank-
Based estimator and the MM-Estimator. This effect is even more pronounced when the
distribution has heavy tails.

24

Chapter 4

Comparing Regression Slopes with Robust Methods

4.1 Background

In social sciences and many other studies, a common interest is to determine whether
associations between two variables differ across groups. For example, does the
association between antisocial behavior and cognitive ability differ across races? Does
the association between the risk of lung cancer and smoking differ across gender? Or as
we explored in the applications to the cortisol and wellbeing measures section of this
thesis, does the association between the cortisol awakening response and some key
measures of wellbeing differ before and after an intervention.

25

4.2 Comparing Regression Slopes with Bootstrap Method

To compare the slopes of two groups, this thesis applied the percentile bootstrap
technique to estimate the confidence intervals of the difference between two slopes. To
elaborate, for two independent groups, let 𝑦 𝑖𝑗
, 𝑋 𝑖𝑗
be the ith vector of observations in the
jth group, 𝑖 =1,…,𝑛 𝑗 . Suppose

𝑦 𝑖𝑗
=𝛽 0𝑗 +𝑋 𝑖𝑗
′
𝛽 𝑗 +𝜆 𝑗 (𝑋 𝑗 )𝜀 𝑖𝑗
,
where 𝛽 𝑗 =(𝛽 1𝑗 ,…,𝛽 𝑝𝑗
)
′
is a vector of slope parameters for the jth group, 𝜆 𝑗 (𝑋 𝑗 ) is some
unknown function of 𝑋 𝑗 , and 𝜀 𝑖𝑗
has variance 𝜎 𝑗 2
. The problem considered here is to
compute a 1−𝛼 confidence interval for 𝛽 𝑘 1
−𝛽 𝑘 2
, 𝑘 =1,…,𝑝 , the difference between
the slope parameters associated with the kth predictor. It is not assumed that 𝜖 𝑖 1
and 𝜖 𝑖 2

have a common variance. Moreover, the goal is to get an accurate confidence interval
without specifying what λ might be.
The conventional method for comparing slope parameters (e.g., Huitema, 1980)
performs poorly when standard assumptions are violated. Conerly and Mansfield (1988)
provide references to other solutions. Included is Chow's (1960) likelihood ratio test,
which is also known to fail. (For results related to the method described here, see Wilcox
1997c.)
When using any robust estimator with a reasonably high breakdown point, the
percentile bootstrap technique appears to give reasonably accurate confidence intervals

26

for a fairly broad range of nonnormal distributions and heteroscedastic error terms. This
suggests a method for addressing the goal considered here, and simulations support its
use.
Briefly, the procedure begins by generating a bootstrap sample from the jth group.
For the jth group, randomly sample 𝑛 𝑗 vectors of observations, with replacement, from
(𝑦 1𝑗 ,𝑋 1𝑗 ),…,(𝑦 𝑛 𝑗 𝑗 ,𝑋 𝑛 𝑗 𝑗 ) . Let 𝑑 𝑘 ∗
=𝛽 ̂
𝑘 1
∗
−𝛽 ̂
𝑘 2
∗
be the difference between the resulting
estimates of the kth predictor, 𝑘 =1,…,𝑝 . Repeat this process B times yielding
𝑑 𝑘 1
∗
,…,𝑑 𝑘𝐵
∗
. Put these B values in ascending order yielding 𝑑 𝑘 (1)
∗
≤⋯≤𝑑 𝑘 (𝐵 )
∗
. Let 𝑙 =
𝛼𝐵 /2, 𝑢 =(1−
𝛼 2
)𝐵 , rounded to the nearest integer, in which case an approximate 1−
𝛼 confidence interval for 𝛽 𝑘 1
−𝛽 𝑘 2
is (𝑑 𝑘 (𝑙 +1)
∗
,𝑑 𝑘 (𝑢 )
∗
) .
The R function reg2ci computes a 0.95 confidence interval for the difference
between regression slope parameters corresponding to two independent groups using the
bootstrap method just described. The default number of bootstrap samples is 599, which
is what was used in this simulation study.

4.3 Simulation Design – Test of the Robust Methods when Comparing
Regression Slopes

Our simulation experiment consists of two parts. In part 1, we focus on testing the
estimators' performances in controlling type I error probabilities. Normal, skewed, heavy-

27

tailed and both skewed and heavy-tailed distributions in independent variable X and error
term 𝜖 were simulated. Both homoscedasticity and heteroscedasticity were simulated as
well.
In part 2, we look at the statistical power of those estimators under different
simulated distribution scenario – four distributions (normal, skewed light-tailed,
symmetric heavy-tailed and skewed heavy-tailed) were simulated for X and 𝜖 and two
variance patterns (homoscedasticity and heteroscedasticity) were simulated as well.

4.3.1 Simulation Part 1 – Test of Type I Error Probability

The goal of this part of simulation is to identify the impact of non-normality and
heteroscedasticity on the performances of the different methods in terms of Type I error
control. Data are generated from the model

𝑌 𝑗 =𝑋 𝑗 +𝜆 (𝑋 𝑗 )𝜀 𝑗 , 𝑗 =1,2 ,

where λ is a function of 𝑋 𝑗 used to model heteroscedasticity. We focus on the situation of
two independent groups. The independent variable X and the error term 𝜖 are generated
from the g-and-h distributions (Hoaglin, 1985), which contain the standard normal
distribution as a special case. If Z has a standard normal distribution, then

28

𝑊 =
{

exp(𝑔𝑍 )−1
𝑔 exp(
ℎ𝑍 2
2
), 𝑖𝑓 𝑔 >0
𝑍𝑒𝑥𝑝 (
ℎ𝑍 2
2
), 𝑖𝑓 𝑔 =0

has a g-and-h distribution.
When 𝑔 =0 and ℎ=0, W has a standard normal distribution. Skewness and
heavy-tailedness of the g-and-h distribution are determined by the values of g and h,
respectively. As the value of g increases, the distribution becomes more skewed. As the
value of h increases, the distribution becomes more heavy-tailed.
Four types of distributions are considered for 𝑋 𝑗 : standard normal (g=0, h=0),
asymmetric light-tailed (g=0.2, h=0), symmetric heavy-tailed (g=0, h=0.2), and
asymmetric heavy-tailed (g=0.2, h=0.2). More properties of the g-and-h distribution can
be found in Table 4-1. The error term 𝜖 𝑗 is also randomly generated based on one of these
four g-and-h distributions. Additional properties of the g-and-h distribution are
summarized by Hoaglin (1985). To simulate heteroscedasticity, we let 𝜆 (𝑋 )=|𝑋 |+1.

29

Table 4-1. Some properties of the g-and-h distribution
g h 𝜅 1
𝜅 2

0.0 0.0 0.00 3.00
0.0 0.2 0.00 21.46
0.2 0.0 0.61 3.68
0.2 0.2 2.81 155.98

The following sample sizes are considered:

1. 𝑛 1
=20,𝑛 2
=20 ;
2. 𝑛 1
=20,𝑛 2
=40 ;
3. 𝑛 1
=40,𝑛 2
=40 .

We purposely choose to keep 𝑛 1
the same while increasing 𝑛 2
in order to examine
whether the performance of the methods can be improved by increasing the sample size
of one group. Theoretically, the performance of a method improves as the sample sizes of
the groups increase. However, in many situations, it is easier to obtain additional samples
from one population but not the other. We are interested in understanding whether
additional samples in one group could improve the accuracy of a method.

30

Table 4-2 provides a brief summary of the simulation design. A total of 24
simulated conditions are considered and 6 methods are tested, so there are altogether 144
simulation tests. (One note is that the MM-estimation does not converge when sample
sizes are relatively small, so we are only able to obtain the results of MM-estimation
when both sample sizes are 40.) The probability of a Type I error is based on 2,000
replications. The actual Type I error probability, when testing at the 𝛼 =.05 level, is
estimated with 𝛼̂ , the proportion of p values less than or equal to .05.

Table 4-2. Summary of the simulation design
𝑿 and 𝒆 VP (Variance Pattern) Sample Size
Standard normal
Skewed light-tailed
Symmetric heavy-tailed
Skewed heavy-tailed
𝜆 (𝑋 𝑖 )=1
𝜆 (𝑋 𝑖 )=|𝑋 𝑖 |+1
𝑛 1
=20, 𝑛 2
=20
𝑛 1
=20, 𝑛 2
=40
𝑛 1
=40, 𝑛 2
=40

4.3.2 Simulation Part 2 – Power Analysis

In this part of simulation, we compare the methods' statistical power. Data are generated
from the model

31

𝑌 1
=𝑋 1
+𝜆 (𝑋 1
)𝜖 1
,
𝑌 2
=𝜆 (𝑋 2
)𝜖 2
.

Similar to Simulation Part 1, λ is a function of X used to model heteroscedasticity.
Data are generated from the g-and-h distributions. Four types of distributions are
considered for both X and 𝜖 : standard normal (𝑔 =0, ℎ=0), asymmetric light-tailed
(𝑔 =0.2, ℎ=0), symmetric heavy-tailed (𝑔 =0, ℎ=0.2), and asymmetric heavy-tailed
(𝑔 =0.2, ℎ=0.2). To simulate heteroscedasticity, we let 𝜆 (𝑋 )=|𝑋 |+1.

4.4 Simulation Results

4.4.1 Simulation Results of the Actual Type I Error Rates

This section reports the results of the simulation tests in Section 4.3.1. The performances
of the 6 methods considered in this thesis for comparing regression slopes in terms of
control of Type I error probability are reported and compared. Results are shown in Table
4-3 to Table 4-8.

32

4.4.1.1 Type I Error Rates with Small Sample Sizes

As shown in Table 4-3, when both sample sizes are 20 and there is no heteroscedasticity
in the generated data, all methods except LTS (i.e., Rank-Based, Theil-Sen, Modified
Theil-Sen and OLS) provide reasonable control over Type I errors in all four distribution
conditions. (LTS performs poorly across all simulations and we will only show its results
in the tables but will not report its performance verbally in this section.)
More specifically, all methods but LTS yield actual Type I error reasonably close
to the nominal level .05. Actually the average Type I error rates estimated by the four
methods range from .028 to .067, which meets the criterion suggested by Bradley (1978)
that as a general guide, at a minimum the actual Type I error probability should be
between .025 and .075 when testing at the .05 level.

Table 4-3. Actual Type I error rates (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟐𝟎 ; with homoscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.039
.049
.040
.041
.034
.033
.031
.030
.030
.031
.029
.028
.00
.00
.001
.0005
.067
.063
.061
.059

33

When we take a closer look at Table 4-3, we notice that the actual Type I errors
estimated by the OLS method go above the nominal level while the alternative robust
methods give Type I errors that drop below the nominal level. Relatively the Theil-Sen
method and the Modified Theil-Sen method are more conservative whose actual Type I
errors range from .028 to .034. But as we mentioned in the previous paragraphs, all four
methods perform well in the conditions of Table 3-3 by providing the actual Type I error
rates close enough to the nominal level.

Table 4-4. Actual Type I error rates (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟐𝟎 ; with heteroscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.055
.059
.056
.048
.046
.044
.044
.042
.041
.042
.038
.041
.002
.0025
.0015
.001
.084
.090
.087
.091

When there is heteroscedasticity and in the meanwhile we still keep the same
sample sizes as in Table 4-3, that is, 𝑛 1
=𝑛 2
=20, from the results shown in Table 4-4
we can see that the three robust methods (Rank-Based, Theil-Sen and Modified Theil-
Sen) continue to perform well but the OLS method deteriorates in terms of the control
over the probability of Type I errors. In more details, when both sample sizes are 20 and

34

there is heteroscedasticity in the generated data, the actual Type I error rates given by the
three robust methods range from .038 to .059 while the actual Type I error rates given by
the OLS method are between .084 to .091 which deviate from the nominal level.

4.4.1.2 Type I Error Rates with Increased Sample Sizes

If we increase one of the sample sizes to 40, the three robust methods still provide good
control over the Type I error rates under both heteroscedasticity and homoscedastity. As
for the OLS method, additional samples do seem to have an effect on its performance. As
shown in Table 4-6, the actual Type I error rates given by the OLS method range from
.072 to .079, which, compared to Table 4-4, is an improvement. However, the estimated
Type I error rates are still above .075 when there are heavy tails, with actual Type I error
rates around .078.

35

Table 4-5. Actual Type I error rates (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.043
.051
.041
.041
.033
.033
.033
.032
.032
.030
.030
.030
.002
.001
.001
.0005
.064
.061
.062
.058

Table 4-6. Actual Type I error rates (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.048
.058
.056
.055
.043
.040
.037
.035
.042
.041
.039
.040
.0015
.002
.002
.002
.075
.072
.079
.078

36

If we increase both sample sizes to 40, then the OLS method provides good
control over the probability of a Type I error when there is homoscedasticity in the
generated data (as shown in Table 4-7). When there is heteroscedasticity, the
performance of the OLS method continues to improve but still does not deliver
consistently well when there are heavy tails (see Table 4-8).
Table 4-7 and Table 4-8 also report the Type I error rates given by the MM-
estimator. As we can see, MM-estimator performs well under the four simulated
distributions for both the heteroscedasticity and homoscedasticity variance patterns.

Table 4-7. Actual Type I error rates (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.053
.053
.052
.048
.031
.032
.050
.029
.033
.032
.029
.028
.039
.044
.030
.036
.0065
.0045
.0035
.0030
.054
.067
.052
.055

37

Table 4-8. Actual Type I error rates (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.053
.054
.053
.054
.040
.035
.036
.036
.038
.033
.036
.036
.031
.032
.024
.030
.0025
.0060
.0050
.0050
.070
.068
.078
.072

4.4.2 Power Analysis Results

We consider 4 distributions (standard normal, skewed light-tailed, symmetric heavy-
tailed, skewed heavy-tailed), 2 variance patterns - heteroscedasticity and
homoscedasticity, and 3 sets of sample sizes (𝑛 1
=𝑛 2
=20; 𝑛 1
=20 and 𝑛 2
=40; 𝑛 1
=
𝑛 2
=40).

4.4.2.1 Power Analysis Results with Small Sample Sizes

When both sample sizes are 20 and there is homoscedasticity in the generated data (see
Table 4-9), the OLS method has better power than the three robust methods. When there

38

is heteroscedasticity (see Table 4-10), the power given by those methods drop to below
.30.

Table 4-9. Power Analysis Results (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟐𝟎 ; with homoscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.72
.73
.61
.61
.60
.61
.50
.50
.64
.64
.61
.61
.103
.03
.03
.03
.80
.80
.72
.70

Table 4-10. Power Analysis (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟐𝟎 ; with heteroscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.21
.24
.20
.19
.21
.21
.19
.19
.20
.20
.18
.18
.016
.013
.015
.012
.30
.28
.25
.24

39

4.4.2.2 Power Analysis Results with Increased Sample Sizes

If we increase the sample sizes, the power increases as well (see Table 4-11 to Table 4-
14). When both sample sizes are 40 and there is homoscedasticity (see Table 4-13), the
power given by the robust methods is ≥.96 under all four distributions. The OLS method
gives slightly lower power than the robust methods when there are heavy tails.

Table 4-11. Power Analysis Results (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.83
.84
.61
.61
.76
.76
.50
.50
.78
.78
.77
.77
.25
.27
.30
.30
.89
.88
.82
.82

40

Table 4-12. Power Analysis (𝒏 𝟏 =𝟐𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-Based Modified
Theil-Sen
Theil-Sen LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.26
.27
.24
.24
.25
.25
.22
.22
.24
.24
.22
.22
.019
.033
.021
.024
.33
.32
.25
.25

4.4.2.3 Power Analysis Results Under Heteroscedasticity

When there is heteroscedasticity, even if we increase both sample sizes to 40, no method
gives satisfactory power. Actually the power given by the robust methods is below 0.40
under all four distributions and the power given by the OLS method is below 0.45 when
there are no heavy tails and around 0.30 when there are heavy tails (see Table 4-14).

41

Table 4-13. Power Analysis Results (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.98
.98
.97
.97
.96
.96
.96
.96
.96
.96
.96
.96
.97
.97
.97
.97
.62
.64
.70
.71
.99
.98
.94
.93

Table 4-14. Power Analysis (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
.0
.2
0
.2
.0
.0
.2
.2
.37
.36
.31
.32
.38
.39
.36
.36
.38
.38
.35
.36
.29
.30
.25
.25
.076
.077
.068
.085
.44
.43
.31
.30

42

4.5 Discussion

In this chapter, we have examined the performance of the conventional and alternative
tests for comparing regression slopes when the normality and homoscedasticity
assumptions are violated. When both sample sizes are small, i.e., 20, and there is
homoscedasticity in the generated data, all methods except LTS and MM-Estimator, i.e.,
the Rank-Based method, the Theil-Sen method, the Modified Theil-Sen method and the
OLS method, provide reasonable control over Type I errors under all four simulated
distributions. (The LTS method performs poorly across all simulations and the MM-
estimator has a converging issue when sample size is small.) When there is
heteroscedasticity in the generated data, the three robust methods (Rank-Based, Theil-Sen
and Modified Theil-Sen) still perform well but the OLS method deteriorates in terms of
the control over the probability of Type I errors.
Another goal is to explore whether increasing the sample size of only one group
would affect the performance of a method. When we increased one of the sample sizes to
40, the three robust methods still provide good control over the Type I error rates under
both heteroscedasticity and homoscedasticity variance patterns. The additional samples
do seem to have an effect on the OLS method's performance in which it gives reasonable
actual Type I error rates under heteroscedasticity in some cases (normal distribution and
skewed light-tailed distribution). However, when there are heavy tails in the generated

43

data, the Type I error rates given by the OLS method are still unsatisfactory. When we
increased both sample sizes to 40, the performance of the OLS method continues to
improve but still does not perform consistently well when there are heavy tails. As for the
MM-Estimator, when it does converge (in our case, when both sample sizes are 40), it
provides good control over Type I error probability.

44

Chapter 5

ANCOVA: Robustness of Rank-Based and Other Estimators As Well
As Non-Parametric Estimator on Design Points

5.1 Background

The analysis of covariance (ANCOVA) problem is to compare two independent groups
based on some outcome of interest, Y, in a manner that takes into account some covariate,
X. A classic and well-known approach is the OLS method, which assumes that the error
term of the linear regression model is homoscedastic and has a normal distribution, the
regression lines associated with each group are parallel, and the variances associated with
the error terms for each group are assumed to be identical. More formally, if for the jth
group (𝑗 =1,2), then there are 𝑛 𝑗 randomly sampled pairs of observations, say 𝑋 𝑖𝑗
,𝑌 𝑖𝑗
,
𝑖 =1,…,𝑛 𝑗 , the classic assumption is that for the jth group,

𝑌 𝑖𝑗
=𝛽 𝑋 𝑖𝑗
+𝛽 0𝑗 +𝜀 𝑖𝑗
Eq. 5-1

45

where 𝜀 𝑖𝑗
has variance 𝜎 𝑗 2
, 𝜎 1
2
=𝜎 2
2
, and 𝜖 𝑖𝑗
is independent of 𝑋 𝑖𝑗
. So by implication, for
each group, the conditional variance of Y, given X, does not vary with X, and each group
has the same slope.
It is known that violating one or more of these assumptions can result in serious
practical problems. Concerns about the robustness of the method date back to at least
Atiqullah (1964) who concluded that non-normality is a practical problem. Another
obvious concern is the assumption that the regression lines are parallel. There are several
robust methods for testing this assumption (e.g., Wilcox, 2003 & 2005), but it remains
unclear when such tests have enough power to detect situations where having non-
parallel lines is a practical concern. Yet another concern about equation (5-1) is the
assumption that the association between Y and X is linear. Many methods have been
derived to eliminate the assumption that the association is linear, however, some of these
methods require homoscedasticity and for most there are few if any simulation results
that support their use with small to moderate sample sizes. A simple and very flexible
approach to ANCOVA is described in Wilcox (2003, section 14.8). It allows nonlinearity
and heteroscedasticity, and in the event standard assumptions are met, it has nearly the
same amount of power as the classic ANCOVA method (e.g., Wilcox, 2005, p. 526).
Roughly, the method performs multiple comparisons based on robust measures of
location in a manner that controls the probability of one or more Type I errors.

46

5.2 ANCOVA with the Goal of Comparing Groups in Terms of Some
Measure of Location

There is a vast literature on the analysis of covariance where the goal is to compare
groups in terms of some measure of location, usually means, while taking into account
some covariate. For an entire book devoted to the subject, see Huitema (1980). For a
review of some recent developments, see Rutherford (1992) and Harwell (2003). We are
not able to describe all relevant methods here. Rather, attention is paid to situations where
comparisons are to be made using some robust measure of location.
For the jth group, let 𝑚 𝑗 (𝑥 ) be some population measure of location associated
with y given x. Given x, the problem is determining how the typical value of y in the first
group compares to the typical value in the second. For example, the measure could be the
20% trimmed mean.
The most common strategy is to assume that a straight regression line can be used
to predict y given x . That is, assume that for the jth group, 𝑚 𝑗 (𝑥 )=𝛽 1
𝑥 1𝑗 +𝛽 0𝑗 , 𝑗 =
1,2. Next, estimate the slope and intercept for each group, and use the estimates to
compute a confidence interval for 𝑚 1
(𝑥 )−𝑚 2
(𝑥 ) . Then a percentile bootstrap procedure
will provide probability coverage when working with some robust regression method.
Here the plan is to use the five robust regression methods mentioned earlier as well as the

47

traditional OLS regression method to estimate the slope and intercept for the regression
part.

5.3 Methods Based on Specific Design Points

5.3.1 A Non-Parametric Method

What happens if a regression straight line is not assumed? In this section we will describe
a method for computing a 1−𝛼 confidence interval for 𝑚 1
(𝑥 )−𝑚 2
(𝑥 ) that makes no
parametric assumption about how y is related to x. Not only the a regression straight line
is not assumed, but the assumption of complete heteroscedasticity is allowed, meaning
that the error term for each group can be heteroscedastic, and nothing is assumed about
how the variance of the error term in the two groups are related. Also, as we mentioned in
Section 5.2, 𝑚 𝑗 (𝑥 ) could be any measure of location or scale, but for now we focus on
the 20% trimmed mean.
The general strategy is to use a running interval smoother to approximate the
regression lines and then use the components of the smoother to make comparisons at
appropriate design points. To elaborate, first assume that an x is chosen with the goal of
computing a confidence interval for 𝑚 1
(𝑥 )−𝑚 2
(𝑥 ) . For the jth group, let 𝑥 𝑖𝑗
, 𝑖 =

48

1,…,𝑛 𝑗 be values of the predictors that are available. As previously indicated, 𝑚 𝑗 (𝑥 ) is
estimated with the trimmed mean of the 𝑦 𝑖𝑗
values such that i is an element of the set

𝑁 𝑗 (𝑥 )={𝑖:|𝑥 𝑖𝑗
−𝑥 |≤𝑓 𝑗 ×MADN
𝑗 } .

That is, for fixed j, estimate 𝑚 𝑗 (𝑥 ) using the 𝑦 𝑖𝑗
values corresponding to the 𝑥 𝑖𝑗
values
that are close to x. Here 𝑓 𝑗 is the span, that is, the value used by the running interval
smoother. Some studies show that the choice 𝑓 𝑗 =0.8 or 𝑓 𝑗 =1 generally gives good
results.
Let M
j
(x) be the cardinality of the set N
j
(x) . That is, M
j
(x) is the number of
points in the jth group that are close to x, which in turn is the number of y
ij
values used to
estimate m
j
(x) . When m
j
(x) is the 20% trimmed mean of y, given x, the two regression
lines are defined to be comparable at x if M
1
(x)≥12 and M
2
(x)≥12. The idea is that if
the sample sizes used to estimate m
1
(x) and m
2
(x) are sufficiently large, then a
reasonably accurate confidence interval for m
1
(x)−m
2
(x) can be computed. There is a
selection of methods that can be used and among them, Yuen’s method often gives
satisfactory results, and the bootstrap-t is even better, albeit more costly in terms of
computing time.

49

5.3.2 Comparison at Design Points with Robust Regression Estimators

The robust regression estimators that are applied to compare two groups at design points
are the robust methods described in Chapter 2. A bootstrap method with B=100 is used to
estimate the standard error of the estimated dependent values 𝑦̂ at design points. In more
detail, suppose 𝑦̂
1
𝑟 and 𝑦̂
2
𝑟 are respectively the estimated dependent values of group 1 and
group 2 at design points. Let 𝑉𝐴𝑅 (𝑦̂
1
𝑟 ) and 𝑉𝐴𝑅 (𝑦̂
2
𝑟 ) be, respectively, the estimated
variance of the dependent variable of group 1 at design points and the estimated variance
of the dependent variable of group 2 at design points with the bootstrap method. Then
compute the statistic value (𝑦̂
1
𝑟 −𝑦̂
2
𝑟 )/√𝑉𝐴𝑅 (𝑦̂
1
𝑟 )+𝑉𝐴𝑅 (𝑦̂
2
𝑟 ) which is assumed to have
a standard normal distribution. The confidence interval and p-value are obtained
accordingly.
The design points are usually selected according to the researcher’s interest and
the goal of the study. In a simulation test, we can use the method described in Section
5.3.4 to select the design points.

50

5.3.3 Comparison at Multiple Design Points – Hochberg’s Method

When there are multiple design points, it means that we need to make multiple
comparisons. The goal is to test C hypotheses such that the probability of one or more
Type I errors is at most 𝛼 , a simple way of proceeding is to use the Bonferroni method,
meaning that each test is performed at the 𝛼 /𝐶 level. However, several improvements on
the Bonferroni method have been published that are designed to ensure that the
probability of at least one type I error does not exceed some specified value, 𝛼 . One of
the methods is proposed by Hochberg in 1988.
Hochberg’s method for controlling the probability of one or more Type I errors is
applied as follows. Let p
1
,…,p
C
be the p-values associated with the C tests that are
performed, put these p-values in descending order, and label the results p
[1]
≥p
[2]
≥
⋯≥p
[C]
. Beginning with k=1 (step1), reject all null hypotheses if

𝑝 [𝑘 ]
≤𝛼 /𝑘 ,

where 𝛼 is the desired probability of one or more Type I errors. That is, reject all null
hypotheses if the largest p-value is less than or equal to 𝛼 . If 𝑝 [1]
>𝛼 , proceed as
follows:

1. Increment k by 1. If

51

𝑝 [𝑘 ]
≤𝛼 /𝑘 ,
stop and reject all null hypotheses having a p-value less than or equal 𝑝 [𝑘 ]

2. If 𝑝 [𝑘 ]
>𝛼 /𝑘 , repeat step 1.
3. Repeat steps 1 and 2 until a significant result is obtained or all C hypotheses have
been tested.

5.3.4 Choose Design Points in a Simulation

Suppose it is desired to compare two groups at five x values: 𝑧 1
, 𝑧 2,
𝑧 3,
𝑧 4,
and 𝑧 5
. In
practice, an investigator might have some substantive reason for picking certain design
points, but this process is difficult to study via simulations. For illustrative purposes,
suppose the design points are chosen using the following process.
First, for notational convenience, assume that for fixed j, the x
ij
values are in
ascending order. That is, x
1j
≤⋯≤x
n
j
j
. Suppose z
1
is taken to be the smallest x
i1
value
for which the regression lines are comparable. That is, search the first group for the
smallest x
i1
such that M
1
(x
i1
)≥12. If M
2
(x
i1
)≥12, in which case the two regression
lines are comparable at x
i1
, set z
1
=x
i1
. If M
2
(x
i1
)<12, consider the next largest x
i1

value and continue until it is simultaneously true that M
1
(x
i1
)≥12 and M
2
(x
i1
)≥12.
Let i
1
be the value of i. that is, i
1
is the smallest integer such that that M
1
(x
i1
)≥
12 and M
2
(x
i1
)≥12. Similarly, let z
5
be the largest x value in the first group for which
the regression lines are comparable. That is z
5
is the largest x
i1
value such that

52

M
1
(x
i1
)≥12 and M
2
(x
i1
)≥12. Let i
3
=(i
1
+i
5
)/2, i
2
=(i
1
+i
3
)/2 and i
4
=(i
3
+
i
5
)/2. Round i
2
, i
3
, and i
4
down to the nearest integer and set z
2
=x
i
2
1
, z
3
=x
i
3
1
, and
z
4
=x
i
4
1
.
Finally, consider computing confidence intervals for m
1
(z
q
)−m
2
(z
q
) , q=
1,…,5 by applying Yuen’s method. (In 1974, Yuen derived a method for comparing
trimmed means that is designed to allow unequal winsorized variances.) In more detail,
perform Yuen’s test using the y values for which the corresponding x values are close to
the design point z
q
, and control the probability of at least one type I error among the five
tests by using the critical value given by the five-variate Studentized maximum modulus
distribution.

5.4 Simulation Design

Similar to Section 4.3, our simulation experiment consists of two parts. In part 1, we
focus on testing the estimators' performances in estimating type I error probabilities.
Normal, asymmetric light-tailed, symmetric heavy-tailed and asymmetric heavy-tailed
distributions in independent variable X and error term 𝜖 were simulated. Both scenario of
homoscedasticity and heteroscedasticity were simulated as well.
In part 2, we look at the statistical power of those estimators under different
simulated distribution scenario – four distributions (normal, skewed light-tailed,

53

symmetric heavy-tailed and skewed heavy-tailed) were simulated for X and 𝜖 and two
variance patterns (homoscedasticity and heteroscedasticity) were simulated as well.

5.4.1 Simulation Part 1 – Test of Type I Error Probability

The goal of this part of simulation is to figure out the robustness of the regression
methods introduced in Chapter 2 in terms of their performances of controlling Type I
error. Data are generated from the model

𝑌 𝑗 =𝑋 𝑗 +𝜆 (𝑋 𝑗 )𝜀 𝑗 , 𝑗 =1,2 ,

where λ is a function of 𝑋 𝑗 used to model heteroscedasticity. The independent variable X
and the error term 𝜀 are generated from the g-and-h distributions (Hoaglin, 1985), whose
skewness and heavy-tailedness are determined by the values of g and h, respectively. As
the value of g increases, the distribution becomes more skewed. As the value of h
increases, the distribution becomes more heavy-tailed. When both g and h have values of
0, then it has a standard normal distribution.
Four types of distributions are considered for 𝑋 𝑗 and 𝜀 𝑗 : standard normal ((g=0,
h=0), asymmetric light-tailed (g=0.2, h=0), symmetric heavy-tailed (g=0, h=0.2), and
asymmetric heavy-tailed (g=0.2, h=0.2). More properties of the g-and-h distribution can

54

be found in Table 4-1. Additional properties of the g-and-h distribution are summarized
by Hoaglin (1985).
To simulate homoscedasticity and heteroscedasticity, we let 𝜆 (𝑋 )=1 and
𝜆 (𝑋 )=|𝑋 |+1 respectivley.
The following sample sizes are considered:

1. 𝑛 1
=30,𝑛 2
=30 ;
2. 𝑛 1
=30,𝑛 2
=40 ;
3. 𝑛 1
=40,𝑛 2
=40 .

We purposely choose to keep 𝑛 1
the same while increasing 𝑛 2
in order to examine
whether the performance of the methods can be improved by increasing the sample size
of one group.
Also note that we do not consider sample size =20 because situations arise where
five design points cannot always be found for which the regression lines are comparable.
Table 5-1 provides a brief summary of the simulation design. A total of 24
simulated conditions are considered and 7 methods (5 robust methods, the OLS method
and the non-parametric method) are tested, so there are altogether 144 simulation tests.
The probability of a Type I error is based on 2,000 replications. The actual Type I error
probability, when testing at the 𝛼 =.05 level, is estimated with 𝛼̂ , the proportion of p
values less than or equal to .05.

55

Table 5-1. Summary of the simulation design
𝑿 and 𝒆 VP (Variance Pattern) Sample Size
Standard normal
Skewed light-tailed
Symmetric heavy-tailed
Skewed heavy-tailed
𝜆 (𝑋 𝑖 )=1
𝜆 (𝑋 𝑖 )=|𝑋 𝑖 |+1
𝑛 1
=30, 𝑛 2
=30
𝑛 1
=30, 𝑛 2
=40
𝑛 1
=40, 𝑛 2
=40

5.4.2 Simulation Part 2 – Power Analysis

In this part of simulation, we compare the methods' statistical power. Data are generated
from the model

𝑌 1
=𝑋 1
+𝜆 (𝑋 1
)𝜖 1
,
𝑌 2
=𝜆 (𝑋 2
)𝜖 2
.

Similar to Simulation Part 1, H is a function of X and taken to be 𝜆 (𝑋 )=1 and
𝜆 (𝑋 )=|𝑋 |+1 respectively to simulate homoscedasticity and heteroscedasticity. X and
𝜖 are generated from the g-and-h distributions: standard normal (𝑔 =0, ℎ=0),
asymmetric light-tailed (𝑔 =0.2, ℎ=0), symmetric heavy-tailed (𝑔 =0, ℎ=0.2), and
asymmetric heavy-tailed (𝑔 =0.2, ℎ=0.2).

56

5.5 Simulation Results

5.5.1 Simulation Results of the Actual Type I Error Rates

This section reports the results of the simulation tests in Section 5.4.1. The performances
of the 7 methods (5 robust methods, the OLS method and the non-parametric method) for
comparing regression slopes at design points in terms of control of Type I error
probability are reported and compared. Results are shown in Table 5-2 to Table 5-7.
Before we take a look at the results, remind ourselves of the criterion suggested by
Bradley (1978) is that as a general guide, at a minimum the actual Type I error
probability should be between .025 and .075 when testing at the .05 level.

5.5.1.1 Type I Error Rates with Small Sample Sizes

As shown in Table 5-2, when both sample sizes are 30 and there is no heteroscedasticity
in the generated data, the OLS method and the Non-Parametric method provide
reasonable control over Type I errors in all four distribution conditions, though they are

57

both tend to be conservative. That is, their actual Type I error rates fall in the range of
.025 to .075 but below .05, the nominal Type I error probability.
As far as the robust regression methods, the Theil-Sen and Modified Theil-Sen
methods did not perform satisfactorily with their actual Type I errors all below .025. The
Rank-Based method and the MM-Estimator provide actual Type I error rates right above
.025, so in the range required but tend to be conservative when there is no heavy-
tailedness in the distribution. But when there is heavy-tailedness in the distribution, their
actual Type I rates fall below .025.

Table 5-2. Actual Type I error rates (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.025
.026
.021
.020
.022
.020
.018
.018
.016
.015
.013
.012
.027
.028
.018
.020
.005
.004
.003
.003
.036
.034
.024
.025
.029
.032
.032
.028

When there is heteroscedasticity and in the meanwhile we still keep the same
sample sizes as in Table 5-2, that is, 𝑛 1
=𝑛 2
=30, as shown in Table 5-3, the methods
studied give similar results to the scenario of Table 5-2. Firstly, the OLS method and the
Non-Parametric method still have good control over the Type I error probability.
Secondly, the Theil-Sen methods do not perform satisfactorily with its actual Type I error

58

rates under all four distributions below .02. The Rank-Based, the MM-Estimator and the
modified Theil-Sen method perform well when there is no heavy-tailedness in the
generated data though tend to be conservative, however, they do not perform
satisfactorily when there is heavy-tailed ness with their actual Type I error rates falling
below .025.

Table 5-3. Actual Type I error rates (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.028
.029
.019
.019
.025
.024
.017
.018
.015
.018
.012
.011
.033
.034
.024
.022
.007
.008
.008
.009
.052
.052
.045
.044
.030
.038
.030
.030

5.5.1.2 Type I Error Rates with Increased Sample Sizes

To see if additional samples have a notable effect on the Type I error probability
control of the methods studied, we increase one of the sample sizes to 40. The results are
shown in Table 5-4 and Table 5-5. From the results shown, we did not see a remarkable
improve in the Type I error probability control of the robust regression methods since
their actual Type I error rates are similar to the results when the sample sizes are not

59

increased. That is, the actual Type I error rates provided the Theil-Sen and the modified
Theil-Sen methods are below .025 under four distributions, the Rank-Based and the MM-
Estimator do give actual Type I error rates between .025 and .033 when there is no
heavy-tailedness but do not perform well when there is heavy-tailedness in the generated
data, under either homoscedasticity or heteroscedasticity.

Table 5-4. Actual Type I error rates (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.032
.026
.021
.023
.021
.024
.019
.020
.016
.017
.015
.017
.032
.033
.023
.020
.005
.005
.007
.004
.040
.039
.033
.033
.036
.038
.030
.032

60

Table 5-5. Actual Type I error rates (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.035
.027
.022
.023
.026
.028
.023
.023
.023
.023
.017
.020
.036
.033
.026
.022
.008
.006
.007
.007
.053
.050
.048
.044
.037
.035
.031
.028

If we continue to increase the sample size so both sample sizes to be 40 (see Table
5-6 and Table 5-7), we noticed that both the modified Theil-Sen method and the MM-
Estimator provide satisfactory control of Type I error probability under all four
distributions and both homoscedasticity and heteroscedasticity. However, the Rank-
Based method still cannot handle the heavy-tailedness situation well, with the actual
Type I error rates below .025 under homoscedasticity. Interestingly, when there is
heteroscedasticity in the simulated distribution, the Rank-Based and the Theil-Sen
methods perform well, with their actual Type I error rates ranging between .024 to .034.

61

Table 5-6. Actual Type I error rates (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.032
.027
.021
.022
.034
.034
.031
.029
.022
.022
.019
.019
.034
.033
.028
.032
.012
.010
.010
.008
.040
.041
.031
.029
.041
.036
.035
.033

Table 5-7. Actual Type I error rates (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.032
.034
.028
.033
.040
.038
.033
.030
.033
.030
.028
.024
.041
.037
.029
.030
.017
.012
.009
.012
.058
.053
.057
.043
.045
.046
.033
.035

62

5.5.2 Power Analysis Results

We consider 4 distributions (standard normal, skewed light-tailed, symmetric heavy-
tailed, skewed heavy-tailed), 2 variance patterns - heteroscedasticity and
homoscedasticity, and 3 sets of sample sizes (𝑛 1
=𝑛 2
=30; 𝑛 1
=30 and 𝑛 2
=40; 𝑛 1
=
𝑛 2
=40).

5.5.2.1 Power Analysis Results with Small Sample Sizes

When both sample sizes are 30, under homoscedasticity (see Table 5-8), the OLS method
provides the best power and the Non-Parametric method provides the lowest power
across the estimators we studied (except the LTS method, as we mentioned earlier, we
will exclude it in the verbal reports but only show its results in the tables). Among the 4
regression methods, the MM-Estimator gives the best power, especially when there is no
heavy-tailedness in the distribution. Then followed by the Rank-Based method and the
modified Theil-Sen method who perform on par with each other, and then followed by
the Theil-Sen method.

63

Table 5-8. Power Analysis (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.53
.51
.49
.48
.53
.51
.48
.47
.43
.42
.38
.37
.61
.59
.53
.51
.16
.15
.17
.16
.72
.69
.52
.49
.26
.25
.19
.19

When there is heteroscedasticity in the simulated distribution with both sample
sizes still being 30 (see Table 5-9), the OLS method still gives relatively the best power,
but overall the power is low across the methods studied.

Table 5-9. Power Analysis (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.12
.16
.10
.10
.14
.14
.11
.10
.11
.11
.09
.08
.13
.13
.09
.10
.034
.04
.03
.03
.21
.21
.13
.13
.08
.08
.06
.06

64

5.5.2.2 Power Analysis Results with Increased Sample Sizes

If we increase the sample sizes, the power increases as well (see Table 5-10 to 5-13).
When both sample sizes are 40 and there is homoscedasticity in the generated data (see
Table 5-12), the power given by the OLS method is ≥.90 when the distribution is light-
tailed. When there is heavy-tailedness in the data, the MM-Estimator has better power
than the other methods. The Rank-Based method has power lower than the MM-
Estimator but higher than the Theil-Sen method and the modified Theil-Sen method.

Table 5-10. Power Analysis (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.65
.64
.61
.60
.65
.65
.60
.61
.56
.55
.50
.51
.71
.70
.67
.65
.25
.24
.28
.29
.82
.80
.65
.63
.36
.34
.24
.24

65

Table 5-11. Power Analysis (𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.17
.17
.12
.13
.19
.17
.14
.13
.16
.15
.12
.12
.16
.17
.11
.12
.04
.05
.03
.04
.27
.25
.17
.15
.10
.10
.08
.08

5.5.2.3 Power Analysis Results Under Heteroscedasticity

When there is heteroscedasticity in the simulated distribution, even if we increase
both sample sizes to 40, no method gives satisfactory power. The best power we could
get is .30, given by the OLS method under normal distribution. Generally the power
given by the robust regression methods is around .20. The power given by the Non-
Parametric method is the lowest among the methods studied, which is under .15 under all
four distributions simulated.

66

Table 5-12. Power Analysis (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.80
.81
.81
.82
.79
.79
.77
.76
.73
.72
.70
.69
.86
.86
.85
.84
.38
.40
.46
.46
.92
.90
.80
.78
.49
.47
.36
.34

Table 5-13. Power Analysis (𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.22
.21
.16
.17
.23
.23
.20
.20
.21
.21
.18
.17
.21
.21
.16
.15
.056
.067
.056
.058
.30
.29
.20
.19
.14
.13
.10
.09

67

5.6 Concluding Remarks

This chapter studies the robustness of the 5 robust regression methods, the OLS method
and the Non-Parametric method when comparing regression slopes at design points
between two groups. Overall speaking, the OLS method and the Non-Parametric method
perform well under all simulated distributions, however, no robust regression method
performs well in Type I error probability control when sample size is small (n=30). When
sample size is small, that is, 30, the Theil-Sen and modified Theil-Sen method are too
conservative with their actual Type I error rates falling below .025. The Rank-Based
method and the MM-Estimator do provide Type I error rates above .025 and below .05,
however they cannot handle the heavy-tailed distribution well with their Type I error
rates falling below .025. We see similar patterns under both the homoscedasticity and
heteroscedasticity.
When sample size is larger, that is, both sample sizes are 40, the modified Theil-
Sen method and the MM-Estimator provide good control over the Type I error
probability under both homoscedasticity and heteroscedasticity. Interestingly, the Rank-
Based method and the Theil-Sen method perform well under heteroscedasticity but not
so under homoscedasticity when there are heavy tails in the distribution.
As far as the power analysis, the OLS method provides the best power and the
Non-Parametric method provides the worst power across the estimators we studied
(except LTS). However, no method gives good power under heteroscedasticity. More
specifically, even if we increase both sample sizes to 40, the best power we could get is

68

only .30. When there is heavy-tailedness in the distribution, the best power we could get
is only .20. The power given by the Rank-Based method is around .20 when sample sizes
are 40 and under heteroscedasticity.

69

Chapter 6

Impact of Intervention on the Association between Cortisol Awakening
Response and Wellbeing Measures

6.1 Study Background

This study is to determine if an intervention procedure has an effect on the association
between cortisol awakening response (CAR) and three key measures of wellbeing (i.e.,
meaningful activity, perceived control, and life satisfaction). These three key measures
were taken before and after the intervention.
The intervention was a lifestyle-oriented occupational therapy program. Extant
results indicate that the lifestyle intervention program used in the Well Elderly 2 study
tends to improve measures of wellbeing among ethnically diverse elders in community-
based settings, the majority of whom are from populations at risk for health disparities.
Currently, however, little is known regarding what role, if any, cortisol plays in the
effectiveness of this particular intervention strategy in terms of improving psychosocial
measures of wellbeing.

70

6.2 The Participants

The participants were 460 men and women aged 60 to 95 years (mean age 74.9; 66%
female). All participants were residents of, users of, or visitors to the study recruitment
sites, demonstrated no overt signs of psychosis or dementia, and were able to complete
the study assessment battery (with assistance, if necessary). Participants were recruited
from 21sites in the greater Los Angeles area, including nine senior activity centres,
eleven senior housing residences, and one graduated care retirement community. Some
participants (N = 141) did not receive intervention the first six months, but after six
months they received intervention and were measured again six months later. Some
participants (N = 232) did not receive intervention at all before or after the six months.
These two independent groups of participants are the target of this analysis.

6.3 The Intervention

This intervention consisted of small group and individual sessions led by a licensed
occupational therapist. Typically, each group had six to eight members, all recruited from
the same site and treated by the same intervener. Monthly community outings were
scheduled to facilitate direct experience with intervention content such as the use of

71

public transportation. The interveners and participants were blind to the study design and
hypotheses. The intervention has such key elements as: Identification and implementation
of feasible and sustainable activity-relevant changes; Development plans to overcome
mundane obstacles to enacting activity-relevant changes (e.g., bodily aches or
transportation limitations); Participation in selected activities; Rehearsal and repetition of
changes to everyday routine.

6.4 Measures and Variables

6.4.1 Measure of CAR (Cortisol Awakening Response)

Cortisol, the primary hormone product secreted by the hypothalamus-pituitary-adrenal
(HPA) axis, is considered to be a biomarker of HPA axis activity. Dysfunction in the
HPA axis is implicated in the development of a variety of sub-clinical and clinical
conditions.
The change in cortisol upon awakening, and when measured 30-60 min later, is
known as the cortisol awakening response (CAR). One appeal of the CAR is that it is an
easy parameter of HPA axis function to measure because it does not require laboratory
conditions or administration of exogenous agents; rather awakening itself is a consistent,
recurring, and useful index and has since attracted growing attention. Both enhanced and
reduced CARs are associated with various psychosocial factors including measures of

72

wellbeing, depression and anxiety disorders. These results suggest that intervention might
impact the association between the CAR and the psychosocial measures used in the Well
Elderly 2 study. They also suggest that the CAR might play a role in the extent to which
individuals respond to intervention.

6.4.2 Three Key Measures of Wellbeing – MAPA, PC and LSIZ

The frequency of meaningful activities was measured with the Meaningful Activity
Participation Assessment (MAPA) instrument, which was studied by Eakman et al. and
found to be a reliable and valid measure of meaningful activity, incorporating both
subjective and objective indicators of activity engagement. It consists of the sum of 29 7-
point Likert scales.
Perceived Control (PC) was measured with an instrument studied by Eizenman et
al. It consists of the sum of 8 four-point Likert scales. So the possible PC scores range
between 8 and 32. Higher PC scores reflect greater perceived control.
Life satisfaction was assessed by the Life Satisfaction Index-Z (LSIZ), which has
established psychometric properties for use in older adults. The LSIZ consists of 13
items, with responses indicating satisfaction, dissatisfaction, and unsure scored as 2,0,
and 1, respectively. Values on the 13 items are then summed to give a score range of 26
(0 to 26).

73

6.5 Analysis and Results

6.5.1 Analysis of Comparing the Association of CAR and Wellbeing
Before and After the Intervention

The rank-based estimation as well as the other five estimations introduced in the
"Methods" section was applied to the data set and used to estimate the slope of CAR
predicting the well-being measures before and after the intervention. The slopes before
and after the intervention were compared with the bootstrap method introduced in Section
4.2 (loop = 499) and the significance level of whether the slopes before and after the
intervention were different was reported.
Another round of analyses was conducted with outliers among the independent
variable identified and eliminated using the MAD-median rule. The MAD-median rule
for detecting outliers refers to declaring 𝑋 𝑖 an outlier if

|𝑋 𝑖 −𝑀 |
𝑀𝐴𝐷 /.6745
>𝐾 ,

74

where 𝑀𝐴𝐷 is the sample median of the n values |𝑋 1
−𝑀 |,…,|𝑋 𝑛 −𝑀 |, M is the
median of 𝑋 1
,…,𝑋 𝑛 , and K is taken to be √𝜒 0.975,1
2
, the square root of a chi-squared
distribution with one degree of freedom (cf. Davies & Gather, 1993).

6.5.2 Results of Comparing the Association of CAR and Wellbeing
Before and After the Intervention

Table 6-1 to 6-6 report the results of these analyses. The p-values in the tables are the p-
values of comparing the association of the independent variable (CAR) and the dependent
variable before and after the intervention. The row of "Slope 1" report the estimated slope
of the regression of the dependent variable (MAPA, LSIZ and PC respectively) on CAR
before the intervention and "Slope 2" report the estimated slope after the intervention.

75

Figure 6-1. Association between CAR and MAPA before and after the intervention
estimated with the Rank-Based estimator (outliers in IV eliminated)

Table 6-1, 6-3 and 6-5 report the results when the outliers among the independent
variable are not eliminated. In terms of the association between CAR and MAPA, no
methods report that there is a statistically significant difference before and after the
intervention when testing at the .05 level. The Rank-Based, MM-Estimator and the OLS
method captured a significant difference of the association of CAR and LSIZ before and
after the intervention with respectively reported p-values .040, .043 and .027. The Theil-
Sen and the Modified Theil-Sen method report p-values of 0.057 and 0.047 respectively,
which are on the border of the statistical significance level .05. Lastly, if we take a look at
-0.4 -0.2 0.0 0.2 0.4
40 60 80 100 120 140
CAR
MAPA
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Slope
Before Intervention
= 17.01
Slope
After Intervention
= -17.12

Before Intervention
After Intervention +

76

the association between CAR and PC before and after the intervention, no method reports
that there is a statistical significance.

Figure 6-2. Association between CAR and LSIZ before and after the intervention
estimated with the Rank-Based estimator (outliers in IV eliminated)

When the outliers among the independent variable are detected and eliminated
using the MAD-median rule, as shown in Table 6-2, 6-4 and 6-6, the Rank-Based method,
MM-Estimator and the OLS method report a statistical significant difference of the
association of CAR and MAPA before and after the intervention, but no significant
difference is detected before and after the intervention in terms of the association of CAR
and LSIZ or the association of CAR and PC.

-0.4 -0.2 0.0 0.2 0.4
5 10 15 20 25
CAR
LSIZ
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Slope
Before Intervention
= -9.2e-15
Slope
After Intervention
= -6.34
Before Intervention
After Intervention +

77

Figure 6-3. Association between CAR and PC before and after the intervention
estimated with the Rank-Based estimator (outliers in IV eliminated)

Table 6-1. Compare association of CAR and MAPA before and after intervention
(outliers in IV included)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.28
-1.40
-15.56
0.25
0.033
-14.42
0.25
0.00
-14.16
0.20
-0.23
-19.01
0.69
-0.22
-0.81
0.27
-0.31
-12.62
-0.4 -0.2 0.0 0.2 0.4
15 20 25 30
CAR
PC
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+
+
+
+
+
+ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+
+
+
+
+ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+
+
+
+
+
+
+
Slope
Before Intervention
= -9.2e-15
Slope
After Intervention
= -1.3e-13
Before Intervention
After Intervention +

78

Table 6-2. Compare association of CAR and MAPA before and after intervention
(outliers in IV excluded)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.02
17.01
-17.12
0.08
12.70
-15.96
0.07
12.72
-15.69
0.04
16.77
-18.78
0.17
17.84
-11.00
0.04
16.70
-16.56

Table 6-3. Compare association of CAR and LSIZ before and after intervention
(outliers in IV included)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.040
0.00
-4.97
0.047
-3.77
-5.72
0.057
0.00
-5.68
0.043
0.57
-5.05
0.31
-0.31
-4.3
0.027
0.84
-5.4

79

Table 6-4. Compare association of CAR and LSIZ before and after intervention
(outliers in IV excluded)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.20
-9.2e-15
-6.34
0.20
-1.2e-6
-6.71
0.23
0.00
-6.76
0.19
-0.29
-6.08
0.46
-1.45
-6.15
0.11
-0.055
-6.79

Table 6-5. Compare association of CAR and PC before and after intervention
(outliers in IV included)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.42
0.00
0.00
0.35
-3.02
-0.068
0.65
0.00
0.00
0.34
0.44
-0.92
0.54
0.17
-2.74
0.46
0.15
-0.19

80

Table 6-6. Compare association of CAR and PC before and after intervention
(outliers in IV excluded)
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS
p-value
Slope 1
Slope 2
0.93
-9.2e-15
-1.3e-13
0.52
-3.2e-15
-0.18
0.67
0.00
0.00
0.87
-1.11
-1.74
0.55
-0.68
-3.32
0.99
-1.70
-1.65

6.6 Discussion
This chapter applies the robust regression methods described in Chapter 2 and the
traditional OLS method to a real data set for comparing the association of CAR and
wellbeing measures before and after an intervention. An outlier detecting method is
utilized to detect and eliminate the outliers in the independent variable.
From the results, we can see that some methods captured significant difference
between the association before and after the intervention and some not. Also, conducting
the outlier detection method or not also gives us different results. Also we noticed from
the reported results that the estimated slopes differ dramatically across the methods. For
example, if we look at Table 6-3, although all three methods (Rank-Based, MM-
Estimator and OLS) report p-values <.05, however the estimated slopes by these three
methods are very different. The estimated slope 1 is 0.00 by the Rank-Based method,

81

0.57 by the MM-Estimator and 0.84 by the OLS method. This is a real example which
tells us that applying different method could give us very different results. We should
always be cautious about which method to use because different methods could give us
dramatically different results and we need to choose the one that will not mislead us to
the wrong conclusions.
If we take a further look at the simulation studies in this thesis, one limitation in
this study is that it only considers the scenario that X and ε have the same distributions
but ignores the scenario that X and ε have different distributions. In doing so, since the
outliers among X can reduce the standard errors so the OLS method takes advantage of
that and can gain efficiency as a result. To get a sense of that, we conducted additional
simulation tests and included part of the results in Appendix D. The data are generated in
the same way as we did in Chapter 3 to Chapter 5, and the only difference is that to
simulate a different distribution of X and ε , we keep X standard normal distributed while
ε has heavy-tailed or skewed heavy-tailed distribution. Simulations were done to test the
performances of the robust methods in Chapter 2 in terms of their relative efficiency,
comparing regression slopes between two groups, and comparing two groups at design
points.
What we discovered are mainly two things: firstly, when we look at the relative
efficiency of the robust estimators and the traditional OLS method, as we predicted, since
now the outliers among X are not there, the OLS method cannot take advantage of that
any more. In all of the simulated scenario, the relative efficiency of the OLS method

82

versus the robust estimators are >1, which means that the OLS method is not preferred to
any of the robust estimator in any scenario simulated.
Secondly, the relative performance of the robust methods being tested when X and
ε have different distributions has the similar pattern as when X and ε have the same
distributions. Generally speaking, the methods that perform better under HET when X
and ε have the same distributions still perform better under HET when X and ε have
different distributions (more specifically, when X is normally distributed and ε has heavy
tails or both skewness and heavy tails); the methods that cannot handle heavy tails very
well when X and ε have the same distributions still cannot handle heavy tails in e when X
and ε have different distributions.
For future exploration of this project, the scenario that the two groups have
different distributions can be simulated and the performances of the mentioned methods
can be tested and compared. Another note worth addressing is about the way of
simulating the design points. We used the method described in Section 5.3.4 to select the
design points, but it will be interesting to know what will happen if we use other ways to
select the design points, for example, median, quantiles, etc.

83

Bibliography

Atiqullah, M. (1964). The robustness of the covariance analysis of a one-way
classification. Biometrika, 51, 365–372.
Chang, W., McKean, J. W., Naranjo, J. D., and Sheather, S. J. (1999). High breakdown
rank-based regression, Journal of the American Statistical Association, 94, 205-219.
Clark F, Azen SP and Zemke R, et al. (1997). Occupational therapy for independent-
living older adults: A randomized controlled trial. JAMA, 278, 1321-6.
Clark F, Azen SP and Carlson M, et al. (2001). Embedding health-promoting changes
into the daily lives of independent-living older adults: long-term follow-up of
occupational therapy intervention.J Gerontol B Psychol Sci Soc, 56, 60-3.
Croux, C.,Rousseeuw, P. J., \& Hossjer, O.(1994). GeneralizedS-estimators. Journal of
the AmericanStatistical Association, 89, 1271-1281.
Dahlquist, G., & Bjorck, A. (1974). Numerical methods. Englewood Cliffs, NJ: Prentice
Hall.

84

Dietz, E. J. (1987). A comparison of robust estimators in simple linear regression.
Communications in Statistics-Simulation and Computation, 16, 1209-1227.
Dietz, E. J. (1989). Teaching regression in a nonparametric statistics course. American
Statistician, 43, 35-40.
Eakman AM, Carlson ME and Clark, FA. (2010). The meaingful activity participation
assessment: a measure of engagement in personally valued activities. Int J Aging
Hum Dev., 70, 299-317.
Eizenman, D. R., Nesselroade, J. R., Featherman, D. L. and Rowe, J. W. (1997).
Intraindividual variability in perceived control in an older sample: The MacArthur
successful aging studies. Psychology and Aging, 12, 489-502.
Glasgow RE, Lichtenstein E, Marcus AC. (2003). Why dont we see more translation of
health promotion research to practice? Rethinking the efficacy-to-effectiveness
transition.Am J Public Health, 93, 1261-1267.
Golub, G. H., & van Loan, C. F. (1983). Matrix computations. Baltimore, MD: Johns
Hopkins University Press.
Hastie, T. J., & Tibshirani, R. J. (1990). Generalized addictive models. New York, NY:
Chapman and Hall.

85

Hay, J., Labree, L., Luo, R., et al. (2002). Cost-effectiveness of preventive occupational
therapy for independent-living older adults. J Am Geriatr Soc,50, 1381-8.
Hendricks,J. and Hatch, LR. (2009). Theorizing lifestyle: exploring agency and structure
in the life course. In: Bengston VL, Gans D, Putney N, Silverstein M, eds.
Handbook of Theories of Aging. New York, NY: Springer Publishing Company, 435-
54
Hettmansperger, T.P., and McKean, J.W. (1998). Robust Nonparametric Statistical
Methods, London: Arnold.
Hettmansperger, T.P., and McKean, J.W. (2011). Robust Nonparametric Statistical
Methods, 2nd Ed., Boca Raton, FL: Chapman-Hall.
Hoaglin, D. C. (1985). Summarizing shape numerically: The g-and-h distribution. In D.
Hoaglin, F. Mosteller & J. Tukey (Eds.) Exploring Data Tables Trends and Shapes.
New York: Wiley.
Hogg, R.V., McKean, J.W., and Craig, A.T. (2013). Introduction to Mathematical
Statistics, 7th Edition, Boston: Pearson.
Huitema, B.E. (1980). The Analysis of Covariance and Alternatives. New York: Wiley.
Huitema, B.E. (2011). The Analysis of Covariance and Alternatives: Statistical Methods
for Experiments, Quasi-Experiments, and Single-Case Studies, Hoboken, NJ: Wiley.

86

Hussain, S. S., \& Sprent, P. (1983). Non-parametric regression. Journal of the Royal
Statistical Society, 146, 182-191.
Jackson, J., Mandel, D., Blanchard, J.Carlson, M., Cherry, B.,Azen, S.,Chou, C-P.,
Jordan-Marsh, M., Forman, T., White, B., Granger, D., Knight, B. and Clark, F.
(2009). Confronting challenges in intervention research with ethnically diverse older
adults: the USC Well Elderly II trial. Clinical Trials, 6, 90-101.
Jaeckel, L. A. (1972), Estimating regression coefficients by minimizing the dispersion of
the residuals, Annals of Mathematical Statistics, 43, 1449-1458.
Kloke, J.D. and McKean, J.W. (2012), Rfit: Rank-based estimation for linear models, The
R Journal, 4/2, 57-64.
Kloke, J., McKean, J. W., and Rashid, M. (2009), Rank-based estimation and associated
inferences for linear models with cluster correlated errors, Journal of the American
Statistical Association, 104, 384-390.
Koul, H. L., Sievers, G. L., and McKean, J. W. (1987), An estimator of the scale
parameter for the rank analysis of linear models under general score functions,
Scandinavian Journal of Statistics, 14, 131-141.
Low,G. and Molzahn, AE. (2007). Predictors of quality of life in old age: a cross-
validation study. Res Nurs Health, 30, 141-50.

87

Luh, W., \& Guo, J. (2000). Approximate transformation trimmed mean methods to the
test of simple linear regression slope equality. Journal of Applied Statistics, 27, 843-
857.
Maronna, R. A., \& Yohai, V. J. (1993). Bias-robust estimates of regression based on
projections. Annals of Statistics, 21, 965-990.
Maronna, R. A., \& Yohai, V. J. (2010). Correcting MM estimates for "fat" data sets.
Journal of Computational Statistics and Data Analysis, 54, 3168-3173.
Martin, R. A., Yohai, V. J., \& Zamar, R. H. (1989). Asymptotically min-max bias robust
regression. Annals of Statistics, 17, 1608-1630.
McKean, J.W. (2004). Robust analysis of linear models, Statistical Science, 19(4), 562-
570.
McKean, J.W. and Sheather, S. J. (2009), Diagnostic procedures, Wiley Interdisciplinary
Reviews: Computational Statistics, 1(2), 221-233.
McKean, J. W., Sheather, S. J. and Hettmansperger, T. P. (1994), Robust and High
Breakdown Fits of Polynomial Models, Technometrics, 36, 409-415.
McKean, J. W. and Sievers, G. L. (1989), Rank scores suitable for the analysis of linear
models under asymmetric error distributions, Technometrics, 31, 207-218.

88

McKean, J. W. and Vidmar, T. J. (1994), A comparison of two rank-based methods for
the analysis of linear models, The American Statistician, 48, 220-229.
McKean, J. W., Vidmar, T.J., and Sievers, G. L. (1989), A robust two stage multiple
comparison procedure with application to a random drug screen, Biometrics 45,
1281-1297.
McKean, J. W., Naranjo, J. D., and Sheather, S. J. (1996), Diagnostics to detect
differences in robust fits of linear models, Computational Statistics, 11, 223-243.
McKean, J. W., Naranjo, J. D., and Sheather, S. J. (1999), Diagnostics for comparing
robust and least squares fits, Journal of Nonparametric Statistics, 11, 161-188.
Ng, M. and Wilcox, R. (2010). Comparing the regression slopes of independent groups.
British Journal of Mathematical and Statistical Psychology. 63, pp. 319-340
Olejnik, S. F. and Algina, J. (1985), A review of nonparametric alternatives to analysis of
covariance, Evaluation Review, 9, 51-83.
Peri, K., Kerse, N., Robinson, E., et al. (2008). Does functionally based activity make a
difference to health status and mobility? A randomised controlled trial in residential
care facilities (The Promoting Independent Living Study; PILS). Age Aging, 37, 57-
63.
Puri, M.L. and Sen, P.K. (1969), Analysis of covariance based on general rank scores,
Annuals of Mathematical Statistics, 40, 610–618.

89

Quade, D. (1967), Rank analysis of covariance, Journal of the American Statistical
Association, 62, 1187–1200.
Rashid, M., McKean, J. W. and Kloke, J., D. (2012), R Estimates and associated
inferences for mixed models with covariates in a multi-center clinical trial, Statistics
in Biopharmaceutical Research, 4, 37–49.
Rowe, JW., Kahn, RL. (1998). Successful aging. Aging (Milando), 10, 142-4.
Scholz, F. W. (1978). Weighted median regression estimates. Annals of Statistics, 6, 603-
609.
Schrader, R. M. and McKean, J.W. (1977), Robust analysis of variance, Communications
in Statistics, Part A-Theory and Methods, 6, 879-894.
Scott, JC., Conner, DA., Venohr, I., et al. (2004). Effectiveness of a group outpatient visit
model for chronically ill older health maintenance organization members: a 2-year
randomized trial of the cooperative health care clinic. J AM Geriatr Soc, 52, 1463-
70.
Sen, P. K. (1968). Estimate of the regression coefficient based on Kendall's tau. Journal
of the American Statistical Association, 63, 1379-1389.
Sievers, G. L., (1978). Weighted rank statistics for simple linear regression. Journal of
the American Statistical Association, 73, 628-631.

90

Sievers, G. and McKean, J. W. (1986), On the Robust Rank Analysis of Linear Models
with Non-symmetric Error Distributions, Journal of Statistical Inference and
Planning 13, 215-230.
Smith, J., Borchelt, M., Maier, H., et al. (2002). Health and well-being in the young old
and oldest old. J Soc Issues, 58, 715-32.
Stock,W. A. and Okun, M. A. (1982). The construct validity of life satisfaction among
the elderly. J Gerontol, 37(5), 625-627.
Talwar, P. P. (1991). A simulation study of some non-parametric regression estimators.
Computational Statistics & Data Analysis, 15, 309-327.
Terpstra, J. and McKean, J. W. (2005), Rank-Based Analyses of Linear Models using R,
(2005), Journal of Statistical Software, 14, Issue 7.
Theil, H. (1950). A rank-invariant method of linear and polynomial regression analysis.
Indagationes Mathematicae, 12, 85-91.
Wallace, JI., Buchner, DM., Grothaus, L., et al. (1998). Implementation and effectiveness
of a community-based health promotion program for older adults. J Gerontol A Biol
Sci Med Sci, 53,301-6.
Wilcox, R. R. (1996). Confidence intervals for the slope of a regression line when the
error term has nonconstant variance. Computational Statistics & Data Analysis, 22,
89-98

91

Wilcox, R. R. (1998a). A note on the Theil-Sen regression estimator when the regressor
is random and the error term is heteroscedastic. Biometrical Journal, 40, 261-268.
Wilcox, R. R. (1998b). Simulation results on extensions of the Theil-Sen regression
estimator. Communications in Statistics-Simulation and Computation, 27, 1117-
1126.
Wilcox, R. R. (2003). Applying Contemporary Statistical Techniques Testing. San Diego
CA: Academic Press.
Wilcox, R. R. (2005). An affine invariant rank-based method for comparing dependent
groups. British Journal of Mathematical and Statistical Psychology/British
Psychological Society. pp. p. 33--42.
Wilcox, R. R. (2005). An approach to ANCOVA that allows multiple covariates,
nonlinearity and heteroscedasticity. Educational and Psychological Mesasurement.
Vol. 65, pp. 442-450..
Wilcox, R. R. (2007). ANCOVA: A robust omnibus test based on selected design points.
Journal of Modern Applied Statistical Methods/JMASM. Vol. 5, pp. 14-21.
Wilcox, R. R. (2007). Robust ANCOVA: Some small-sample results when there are
multiple groups and multiple covariates. Journal of Applied Statistics. Vol. 34, pp.
353--364.

92

Wilcox, R. R. (2009). Robust ancova using a smoother with bootstrap bagging. British
Journal of Mathematical and Statistical Psychology. Vol. 62, pp. 427--437.
Wilcox, R. R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd
Ed). San Diego CA: Academic Press.
Wilcox, R. R. (2012). Modern Statistics for the Social and Behavioral Sciences: A
Practical Introduction. CRC Press.
Wilcox, R. R. and Clark, F. (in press). Robust regression estimators when there are tied
values. Journal of Modern and Applied Statistical Methods.
Wilcox, R. R., Granger, D. A., Szanton, S. and Clark, F. (in press). In older adults,
intervention impacts the association between the cortisol awakening response and
three measures of wellbeing: meaningful activities, life satisfaction and perceived
control.
Wood, V., Wylie, ML. and Sheafor, B. (1969). An analysis of a short self-report measure
of life satisfaction: correlation with rater judgements. Journal of Gerontology, 24,
465-469.
Yohai, V. J. (1987). High breakdown point and high efficiency robust estimates for
regression. Annals of Statistics, 15, 642-656.

93

Appendix A: R Functions

A.1 R Function reg2ci

The R function reg2ci
reg2ci (x1,y1,x2,y2,regfun=Rfit,nboot=599,alpha=0.05,ploti=T,xout=F)
computes a 0.95 confidence interval for the difference between regression slope
parameters corresponding to two independent groups. The first two arguments (x1 and
y1) contain the data for the first group, and the data for group 2 are contained in x2 and
y2.
The optional argument “regfun” indicates the regression method to be used. As
we did in this thesis, we set “regfun=Rfit” to call the Rank-Based regression method,
“regfun=MMreg” to call the MM-Estimator, “regfun=tshdreg” to call the modified Theil-
Sen estimator, “regfun=tsreg” to call the Theil-Sen estimator, “regfun=ltsreg” to call the
Least Trimmed Squares estimator and “regfun=ols” to call the OLS estimator. If not
specified, tsreg, the Theil-Sen estimator, is used.
The default number of bootstrap samples, nboot, is B=599, and alpha defaults to
0.05 if unspecified. When set “xout=T”, the MAD-Median method will be used to

94

exclude the outliers in the independent variable. When the argument plotit equals T, the
function also creates a scatterplot that includes the regression lines for both groups.

95

A.2 R Function ancts

The R function ancts
ancts(x1,y1,x2,y2,pts=NULL,regfun=Rfit,nboot=100,alpha=0.05,ploti=T,xout=F)
compares the regression lines of two independent groups at specified design points.
Assume the data are in x1, y1, x2, y2, the design points can be specified by setting “pts”.
For example, if pts contains the values, 5, 8, and 12, the function will compare the
regression lines at x=5,8, and 12.
The optional argument “regfun” indicates the regression method to be used. As
we did in this thesis, we set “regfun=Rfit” to call the Rank-Based regression method,
“regfun=MMreg” to call the MM-Estimator, “regfun=tshdreg” to call the modified Theil-
Sen estimator, “regfun=tsreg” to call the Theil-Sen estimator, “regfun=ltsreg” to call the
Least Trimmed Squares estimator and “regfun=ols” to call the OLS estimator. If not
specified, tsreg, the Theil-Sen estimator, is used.
The default number of bootstrap samples, nboot, is B=100, and alpha defaults to
0.05 if unspecified. When set “xout=T”, the MAD-Median method will be used to
exclude the outliers in the independent variable. When the argument plotit equals T, the
function also creates a scatterplot that includes the regression lines for both groups.

96

A.3 R Function ancova

The R function ancova
ancova(x1,y1,x2,y2,fr1=1,fr2=1,tr=0.2,alpha=0.05,ploti=T,xout=F,pts=NA,sm=F,…)
compares trimmed means at design points with the non-parametric method. The data for
group 1 are stored in x1 and y1, and for group 2, they are stored in x2 and y2. The
arguments fr1 and fr2 are the values of the span for the first and second group, used by
the running interval smoother, which default to 1 if unspecified, and “tr” is the amount of
trimming which defaults to 0.2. The default value for alpha, the probability of at least one
Type I error, is .05. If the argument “pts=NA”, five design points are chosen as described
in Section 5.3.3. The results are returned in “ancova$output”.
If values are stored in “pts”, the function compares groups at the values specified.
So if “pts” contains the values, 5, 8, and 12, the function will test H
0
: m
1
(x)=m
2
(x) at
x=5,8,and 12. When “plotit=T”, the function creates a scatter plot and smooth for both
groups by calling a function which creates a separate smooth for each group. If the
argument “xout=T” is used, leverage points are removed when plotting the regression
lines.

97

Appendix B: 20% Trimmed Mean

In Section 5.3.1, when we describe the Non-Parametric method to estimate the slope at
design points, we mentioned that the measure of location we use in the simulation study
for this non-parametric method is 20%-trimmed mean.
Why 20% trimmed mean? When sampling from a normal distribution, little
efficiency is lost when using the 20% trimmed mean, versus the sample mean. However,
for a heavy tailed distribution, the efficiency of the 20% trimmed mean can be much
higher than the mean. But too much trimming (the maximum amount of trimming is 50%
trimming, i.e., the median) can lead to poor efficiency when sampling from a normal
distribution. The choice of 20% trimming is a compromise between no trimming and the
median (Rosenberger & Gasko, 1983; Wilcox, 1997; Wilcox, 2001).
Another consideration in choosing how much trimming is in terms of hypothesis
testing. Too little trimming results in poor control over the probability of a Type I error,
poor probability coverage of the confidence interval, undesirable power properties, and
relatively low power under departures from normality. However, when sampling from a
normal distribution, too much trimming results in low power. Again, the 20% trimming is
a reasonably good compromise between the mean and the median (Staudte &
Sheather,1990; Wilcox, 2001).

98

Appendix C. True Power for Comparing Groups at Design Points with
Parallel Slopes

We did simulation tests to study the true power of the robust methods and the OLS
method when comparing groups at design points with the generated data having parallel
slopes. The true power test is to reject only when all five tests reject instead of at least
one of the five tests rejects. The data are generated from the model

𝑌 1
=𝑋 1
+𝜆 (𝑋 1
)𝜖 1
,
𝑌 2
=1+𝑋 1
+𝜆 (𝑋 2
)𝜖 2
.

The results are shown in Table D-1 to Table D-6. We can see that the MM-Estimator has
the best power among the robust methods under homoscedasticity. The Rank-Based
method and the Modified Theil-Sen method have similar power when there are no heavy
tails in the distribution, while when there are heavy tails, the Rank-Based method has
better power than the Modified Theil-Sen method, and the Modified Theil-Sen method
has better power than the Theil-Sen method. The OLS method has better power than the
other robust methods but this advantage only exists when the distribution has no heavy
tails. Lastly, when there is heteroscedasticity in the generated data, the power to detect

99

the difference at all five design points is very low across all methods being tested. Lastly,
increasing sample size increases power as well.

Table C-1. True Power (with parallel slopes; 𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0 .0 .44 .44 .32 .53 .10 .65 .28
.2 .0 .43 .42 .33 .53 .12 .64 .29
.0 .2 .45 .39 .30 .51 .15 .48 .22
.2 .2 .44 .39 .30 .51 .15 .45 .22

Table C-2. True Power (with parallel slopes; 𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟑𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0 .0 .058 .068 .059 .045 .006 .10 .09
.2 .0 .058 .065 .060 .053 .008 .10 .10
.0 .2 .041 .055 .050 .033 .009 .061 .087
.2 .2 .046 .052 .050 .037 .009 .056 .088

100

Table C-3. True Power (with parallel slopes; 𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.45
.46
.49
.49
.41
.42
.38
.39
.31
.33
.30
.30
.54
.53
.52
.52
.11
.14
.17
.18
.66
.64
.47
.47
.31
.33
.23
.24

Table C-4. True Power (with parallel slopes; 𝒏 𝟏 =𝟑𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0
.2
.0
.2
.0
.0
.2
.2
.056
.059
.034
.033
.060
.062
.056
.058
.056
.056
.052
.054
.030
.036
.027
.029
.007
.008
.008
.009
.078
.086
.046
.048
.086
.090
.080
.079

101

Table C-5. True Power (with parallel slopes; 𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with
homoscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0 .0 .43 .39 .31 .52 .11 .62 .30
.2 .0 .48 .40 .32 .53 .12 .60 .30
.0 .2 .51 .39 .31 .55 .18 .46 .21
.2 .2 .52 .40 .31 .55 .19 .44 .21

Table C-6. True Power (with parallel slopes; 𝒏 𝟏 =𝟒𝟎 ,𝒏 𝟐 =𝟒𝟎 ; with
heteroscedasticity)
g h Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS OLS Non-
Parametric
.0 .0 .028 .031 .027 .027 .0005 .045 .052
.2 .0 .026 .036 .036 .026 .002 .047 .059
.0 .2 .024 .028 .025 .024 .002 .024 .051
.2 .2 .019 .033 .027 .021 .0015 .022 .048

102

Appendix D. Simulation Results When X and ε Have Different
Distributions

D.1 Relative Efficiency
The simulation results of testing the relative efficiency of the robust methods described in
Chapter 2 are shown in TableD-1 to Table D-6. Here X has standard normal distribution; g 1
and h 1 are used to simulate the distribution of the error term ε. The simulation studies conducted
here have ε either heavy-tailed or skewed heavy-tailed.

Table D-1. Relative Efficiency – RE = OLS/Robust Methods (n=20; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.26
1.19
1.22
1.14
1.22
1.41
1.27
1.18
1.04
0.96

103

Table D-2. Relative Efficiency – RE = OLS/Robust Methods (n=20; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.28
1.26
1.35
1.33
1.34
1.32
1.26
1.25
1.05
1.01

Table D-3. Relative Efficiency – RE = OLS/Robust Methods (n=30; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.28
1.21
1.26
1.18
1.25
1.18
1.29
1.21
1.11
1.03

Table D-4. Relative Efficiency – RE = OLS/Robust Methods (n=30; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.39
1.30
1.48
1.38
1.47
1.38
1.37
1.27
1.14
1.06

104

Table D-5. Relative Efficiency – RE = OLS/Robust Methods (n=40; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.28
1.22
1.25
1.18
1.24
1.18
1.28
1.21
1.13
1.05

Table D-6. Relative Efficiency – RE = OLS/Robust Methods (n=40; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM
Estimator
LTS
.0
.2
.2
.2
1.38
1.31
1.46
1.38
1.46
1.38
1.36
1.28
1.12
1.05

105

D.2 Compare Slopes between Two Groups
The Type I error rates of comparing regression slopes between two groups with the
percentile bootstrap method described in Section 4.2 using the robust methods described
in Chapter 2 and the OLS method are shown in Table D-7 to Table D-12. X has standard
normal distribution, while ε has either heavy-tailed distribution or skewed heavy-tailed
distribution generated from the g-and-h distribution.
Table D-7. Actual Type I Error Rates (𝒏 𝟏 =𝟐𝟎 , 𝒏 𝟐 =𝟐𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
LTS OLS
.0
.2
.2
.2
.040
.043
.033
.031
.031
.029
.001
.001
.066
.073

Table D-8. Actual Type I Error Rates (𝒏 𝟏 =𝟐𝟎 , 𝒏 𝟐 =𝟐𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
LTS OLS
.0
.2
.2
.2
.053
.054
.046
.040
.036
.039
.002
.002
.088
.085

106

Table D-9. Actual Type I Error Rates (𝒏 𝟏 =𝟐𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
LTS OLS
.0
.2
.2
.2
.035
.036
.030
.032
.035
.031
.0005
.001
.063
.066

Table D-10. Actual Type I Error Rates (𝒏 𝟏 =𝟐𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
LTS OLS
.0
.2
.2
.2
.038
.043
.038
.039
.035
.039
.002
.002
.077
.075

Table D-11. Actual Type I Error Rates (𝒏 𝟏 =𝟒𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS
.0
.2
.2
.2
.052
.051
.051
.053
.051
.032
.035
.037
.010
.008
.055
.054

107

Table D-12. Actual Type I Error Rates (𝒏 𝟏 =𝟒𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS
.0
.2
.2
.2
.054
.054
.052
.035
.050
.034
.033
.031
.007
.006
.074
.072

108

D.3 Compare Two Groups at Design Points
The Type I error rates of comparing two groups at design points with the bootstrap
method described in Section 5.3.2 using the robust methods described in Chapter 2, the
OLS method and the Non-Parametric method described in Section 5.3.1 are shown in
Table D-13 to D-18. X has standard normal distribution, while ε has either heavy-tailed
distribution or skewed heavy-tailed distribution generated from the g-and-h distribution.
Table D-13. Actual Type I Error Rates (𝒏 𝟏 =𝟑𝟎 , 𝒏 𝟐 =𝟑𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.019
.017
.016
.017
.012
.013
.024
.023
.004
.004
.028
.030
.027
.022

Table D-14. Actual Type I Error Rates (𝒏 𝟏 =𝟑𝟎 , 𝒏 𝟐 =𝟑𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.019
.022
.017
.017
.011
.012
.028
.027
.007
.010
.040
.038
.026
.022

109

Table D-15. Actual Type I Error Rates (𝒏 𝟏 =𝟑𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.025
.019
.020
.020
.014
.013
.026
.026
.007
.006
.036
.033
.031
.032

Table D-16. Actual Type I Error Rates (𝒏 𝟏 =𝟑𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.023
.025
.023
.023
.016
.019
.031
.027
.007
.007
.040
.039
.030
.031

Table D-17. Actual Type I Error Rates (𝒏 𝟏 =𝟒𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HOM)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.025
.019
.031
.031
.023
.023
.030
.029
.010
.010
.036
.035
.039
.037

110

Table D-18. Actual Type I Error Rates (𝒏 𝟏 =𝟒𝟎 , 𝒏 𝟐 =𝟒𝟎 ; HET)
g
1
h
1
Rank-
Based
Modified
Theil-Sen
Theil-
Sen
MM-
Estimator
LTS OLS Non-
Parametric
.0
.2
.2
.2
.034
.032
.033
.033
.028
.026
.030
.031
.013
.013
.046
.047
.039
.041

Abstract (if available)

Abstract This thesis studied robust regression methods, especially rank-based method and its performance compared to other robust regression methods (the Theil-Sen method, the modified Theil-Sen method, the Least Trimmed Squares method and the MM-Estimator) and the traditional OLS method. Three simulation studies were conducted to determine the efficiency of the methods when estimating the regression slopes, the robustness of the methods when comparing slopes between two groups and the robustness of the methods when comparing slopes at design points between two groups. ❧ Generally speaking, in terms of efficiency, the rank-based method and the MM-Estimator have better efficiency than the other robust regression methods under homoscedasticity. But the Theil-Sen and the modified Theil-Sen method have better efficiency under heteroscedasticity. ❧ As far as comparing regression lines, the results indicate that three robust methods (the Rank-Based method, the Theil-Sen method and the modified Theil-Sen method) perform well in the simulation studies, even when there is heteroscedasticity. Another robust method (the MM-Estimator) performs satisfactorily as well when sample size is relatively large. These four methods studied here are recommended for general use. ❧ When it comes to comparing the slopes at design points, no robust regression method performs well in Type I error probability control when the sample size is relatively small. When the sample size is relatively large, the Theil-Sen and the modified Theil-Sen method perform well but tend to be conservative. Non-Parametric method provides good control over Type I error probability across the distributions and variance patterns simulated

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Comparing dependent groups with missing values: an approach based on a robust method

PDF

Comparing robustness to outliers and model misspecification between robust Poisson and log-binomial models

PDF

Using classification and regression trees (CART) and random forests to address missing data

PDF

Sampling strategies based on existing information in nested case control studies

PDF

Disease risk estimation from case-control studies with sampling

PDF

On the latent change score model in small samples

PDF

Evaluating standard error estimators for multilevel models on small samples with heteroscedasticity and unbalanced cluster sizes

PDF

Non-parametric multivariate regression hypothesis testing

PDF

Exploitation of sparse and low-rank structures for tracking and channel estimation

PDF

Comparing skipped correlations: the overlapping case

PDF

Robust estimation of high dimensional parameters

PDF

Estimation of nonlinear mixed effects mixture models with individually varying measurement occasions

PDF

On the theory and applications of structured bilinear inverse problems to sparse blind deconvolution, active target localization, and delay-Doppler estimation

PDF

Making appropriate decisions for nonnormal data: when skewness and kurtosis matter for the nominal response model

PDF

Outlier-robustness in adaptations to the lasso

PDF

Biometric models of psychopathic traits in adolescence: a comparison of item-level and sum-score approaches

PDF

Two essays on financial econometrics

PDF

Exploring robust aternatives to Pearson's r through secondary analysis of published behavioral science data

PDF

A Bayesian region of measurement equivalence (ROME) framework for establishing measurement invariance

PDF

Generalized linear discriminant analysis for high-dimensional genomic data with external information

Asset Metadata

Creator Ma, Jinxia (author)

Core Title Robustness of rank-based and other robust estimators when comparing groups

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Psychology

Publication Date 04/21/2016

Defense Date 11/13/2015

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag OAI-PMH Harvest,rank-based,robust estimators

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Wilcox, Rand (committee chair), Eckel, Sandrah (committee member), John, Richard S. (committee member), Walsh, David A. (committee member)

Creator Email jinxia.ma@gmail.com,jinxiama@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-238643

Unique identifier UC11278516

Identifier etd-MaJinxia-4336.pdf (filename),usctheses-c40-238643 (legacy record id)

Legacy Identifier etd-MaJinxia-4336.pdf

Dmrecord 238643

Document Type Dissertation

Format application/pdf (imt)

Rights Ma, Jinxia

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA