%% This is simpleRe.tex.
%% H. Seltman and Joel Greenhouse, Summer 2015

% set documentation info, app title, and colors)
\def\author{Joel Greenhouse and Howard Seltman}
\def\title{Review of Simple Regression}
\def\titleColor{Re}
\def\titleFontWeight{bold}
\def\titleFontPercent{125}
\def\titleFont{serif}
\def\headerColor{LightYellow}
\def\backgroundColor{Khaki}
\def\verbatimColor{LightYellow}



\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%% Am I ready tab %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{shinyTab}{Am I ready?}

\h3{Am I ready?}

This app is designed to review simple linear regression for those who have studied it before, but may have forgotten
some of the details.

Please take the pre-test first.  If you cannot answer most of these questions correctly, you should review
the corresponding material from a first course in Statistics before continuing with this app.

\h4{Pre-test}
\shinyQuiz{bank=readyQuiz; addNumbers=TRUE}

\shinyGotoTab{name=Context and Usage; text=Go to the Context and Usage Tab}
\end{shinyTab}


%%%%%%%%%%%%%%%%%%%%%%%%%%%% Context and Usage tab %%%%%%%%%%%%%%%%%%%%%%%%
\begin{shinyTab}{Context and Usage}

\h3{Context and Usage}

Simple linear regression is used in the context of
\begin{enumerate}
\item A single quantitative dependent variable (DV)
\item A single quantitative independent variable (IV)
\item Absence of serial correlation (e.g., a financial time series) and absence of grouping (e.g., students nested within classrooms)
\item A linear relationship between the IV and the DV (or transformations of these), at least within the range of
      IV values studied.
\end{enumerate}

It can be used to analyze samples from both randomized experiments and observational studies, but establishes causality
only the in the former.  Possible goals include learning about the relationship between the IV and the DV and/or predicting
future DV values from known IV values.

\begin{shinyCode}
  NG = 30
  fertilizer = round(rnorm(NG, 100, 30), 2)
  heightGain = round(rnorm(NG, 50+0.3*(fertilizer-100), 1.5), 1)
  growth = data.frame(id=1:NG, fertilizer, heightGain)
\end{shinyCode}

Consider a quick artificial example.  Imagine an experiment where plants are randomized to
have differing amounts of fertilizer (in gms) added, and the growth of the plant (in mm) is
measured after one month.  Here are the first and last few lines of data collected (using the R statistical package):

\begin{shinyCodeOutput}
  print(growth[c(1:3,28:30),], row.names=FALSE)
\end{shinyCodeOutput}

\smallskip
Here is a scatterplot with the "best-fit line":
\begin{shinyPlot}
with(growth, plot(fertilizer, heightGain))
abline(lm(heightGain~fertilizer, data=growth))
\end{shinyPlot}

\smallskip

And here are the simple linear regression results:
\begin{shinyCodeOutput}
  summary(lm(heightGain~fertilizer, data=growth))
\end{shinyCodeOutput}

\smallskip
The rest of this app will show you how to setup and analyze your data and interpret the results.

\shinyGotoTab{name=The Model; text=Jump to the Model tab}

\end{shinyTab}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The Model tab %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{shinyTab}{The Model}
\h3{The Model and EDA}

\h4{Least Squares Regression: Model Formulation and Assumptions}

The simple linear regression \textbf{model} specifies that there is a linear
relationship between $Y$ (the dependent or outcome variable) 
and $X$ (the independent or explanatory variable).  

One way to write the model is:
$$ Y_i = \beta_0 + \beta_1 x_i + \epsilon_i $$
$$ \epsilon_i \sim \mathrm{ Normal}(\mathrm{mean}=0, \mathrm{var}=\sigma^2) $$
$$ \epsilon_i \mathrm{\ independent\ of\ } \epsilon_j$$

This describes the relationship between the DV and IV for each particular subject (the $i^{th}$)
in the population.  Here $\beta_0$ is the population intercept, $\beta_1$ is the 
population slope, and $\epsilon_i$ is the \textbf{error} for the $i^{th}$ subject.  In this
context "error" means unexplained variation (beyond what $x_i$ explains), not "mistakes".
The $\sim$ symbol means "is distributed as".  $\sigma^2$ is called the error variance.

If we apply the tools of mathematical statistics to these statements, we obtain:
$$ \mathrm{E}(Y|x) = \beta_0 + \beta_1 x, \ Y|x \sim \mathrm{Normal}(0,\sigma^2)$$

Here $\mathrm{E}(Y|x)$ means the expected value of $Y$ given a particular
value for $x$, and $Y|x \sim$ means the distribution of $Y$ at a particular
value of $x$.  This second form helps us to 
understand the model and what we expect to see in the data more clearly.  At each fixed value
of $x$ we expect to see a variety of $y$ values, but the mean (expected) value of all of those
$y$ values is $\beta_0 + \beta_1 x$.  And the distribution of the $y$ values around
that mean follows a Gaussian distribution with a constant (not varying with $x$) variance.

Another important reason to look at the model equations is as an aid in understanding
the \textbf{meanings of the parameters}.  By thinking about $x=0$, we can see that
the population intercept, $\beta_0$, is interpreted as the true \textbfem{mean} of the DV when the IV is 
set to \textbfem{zero}.  The interpretation of the population slope, $\beta_1$, is the true
change in the \textbfem{mean} of the DV for any \textbfem{one unit} increase in the IV.  And
the meaning of the error standard deviation, $\sigma$, is the true standard deviation of 
the DV at \textbfem{each specific} value of the IV.

Using either form we can summarize the \textbf{model assumptions} as:

\begin{enumerate}
\item $Y$ is normally distributed \textem{at each fixed x value}.
\item On the average, there is a linear relationship between $Y$ and $X$.
\item The variance of $Y$ is constant and does not change as $X$ changes.
\item The errors are independent of each other.
\end{enumerate}

(Although it may not be clear from the above equations, there is also an important
assumption that the IV is measured without error.  This is often called the
fixed-x assumption.)

If any of these are clearly untenable, we cannot apply simple linear regression.
If the assumptions seem reasonable, then we can proceed with our analysis,
but we will check the assumptions along the way, and possibly conclude that
we need to switch to a different approach.

\medskip
\h4{Exploratory Data Analysis}

When working with any data set, once we have recorded the data, we need to 
perform EDA (Exploratory Data Analysis) before proceeding to a more formal
analysis such as simple linear regression.

The \textbf{goals of EDA} include:
\begin{itemize}
\item checking for data entry errors and outliers
\item becoming more familiar with the individual variables
\item exploring the relationship between the IV and the DV(s)
\item preliminary assumption checking and corresponding model choice
\end{itemize}

We can \textbf{classify EDA} into four types based on whether we are looking at variable individually or in pairs and on
whether we are producing graphs or not:
\begin{enumerate}
\item univariate non-graphical
\item univariate graphical
\item bivariate non-graphical
\item bivariate graphical
\end{enumerate}

In the setting of one quantitative DV and one quantitative IV, \textbf{typical EDA} includes the following.  The examples
are from the above plant growth example.
\begin{enumerate}

\item univariate non-graphical: computation of statistics such as mean, minimum, maximum, and standard deviation of each variable
\begin{shinyCodeOutput}
cat("Summary of 'fertilizer'\n")
print(summary(growth$fertilizer))
cat("\nSummary of 'heightGain'\n")
print(summary(growth$heightGain))
\end{shinyCodeOutput}
\smallskip

\item univariate graphical: a histogram or boxplot or stem-and-leaf plot of each variable
\begin{shinyPlot}
par(mfrow=c(1,2))
hist(growth$fertilizer, xlab="fertilizer (gm)", main="")
hist(growth$heightGain, xlab="height gain (mm)", main="")
\end{shinyPlot}
\smallskip

\item bivariate non-graphical: computation of the correlation coefficient between the two variables
\begin{shinyCodeOutput}
cat("Correlation of 'fertilizer' and 'heightGain'\n")
print(with(growth, cor(fertilizer, heightGain)))
\end{shinyCodeOutput}
\smallskip

\item bivariate graphical: a scatterplot with the DV on the y-axis and the IV on the x-axis
\begin{shinyPlot}
with(growth, plot(heightGain~fertilizer))
\end{shinyPlot}
\end{enumerate}

\smallskip
\textbf{Interpretations:}

The univariate EDA helps to identify outliers and is useful for familiarizing yourself with the data.
If the histogram of the DV shows very few different values, then the Normality assumption is untenable.
If the histogram of the IV shows very few different values, then it will be hard to assess the 
linearity assumption, and another analysis method, e.g., ANOVA, may be more appropriate.

The non-graphical bivariate EDA (correlation statistic) gives a general idea of how tightly the points cluster around a
straight line and the sign shows whether an increase in the IV is associated with a positive or a negative change in the
DV.

The graphical bivariate EDA (scatter plot) helps test the equal variance and linearity assumptions.
If the scatter plot shows a marked difference in the spread of $Y$ at different $X$ values, then 
the equal variance assumption is untenable.  If the general pattern of the relationship between
the DV and the IV is not linear, then the linearity assumption is untenable.

Note that nothing in the EDA helps to check the independent errors assumption.



\shinyGotoTab{name=Estimation and residual analysis; text=Go to Estimation tab}

\end{shinyTab}


%%%%%%%%%%%%%%%%%%%%%%%%%%%% Estimation and residual analysis tab %%%%%%%%%%%%%%%%%%%%%%%%
\begin{shinyTab}{Estimation and residual analysis}

We use the method of \textbf{least
squares} to find estimates of the regression coefficients in the
model: 
\begin{itemize}
\item $\hat{\beta}_0$ is the estimate of the $y$-intercept $\beta_0$, and is an estimate of the average value of $Y$ when $X=0$
\item $\hat{\beta}_1$ is the sample estimate of the population slope, $\beta_1$, and is an estimate of how much the $Y$
variable changes with a 1 unit change in the $X$ variable  
\end{itemize}

We write the least squares line also called the \textbf{fitted} line (or the {\em
fitted values}) as

$$ \hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \cdot x $$

The difference between an observed value $y$ and the fitted value
$\hat{y}$ is called a \textbf{residual}.  It represents how far off the
fitted value ($\hat{y}$) is from the observed value ($y$).  It is also
a sample estimate of the true error, $\epsilon$, for each measurement.  Formally,
the residual (denoted by $e$) is $e = y - \hat{y}$.

In order for inferences about the linear regression coefficients
(i.e., $\beta_0$ and $\beta_1$) to be valid, we need to make certain
assumptions about the probability distribution of $Y$, the response
variable in a linear regression model. These assumptions, as shown in the
previous tab, are normality, equal variance, linearity, and independent errors.

After fitting a least squares model to the data, and computing estimates of 
$\hat{\beta_0}$, $\hat{\beta_1}$, the individual predictions, $\hat{y_i}$,
and the residuals, $e_i$ (or $\hat{\epsilon_i}$), we will use informal
diagnostic plots to provide information on whether the assumptions
of the model are approximately true.  These plots (and the assumptions they check) include
\begin{enumerate}
\item Scatterplot of the residuals ($e$) versus $\hat{y}$ ({\em
linearity, constant variance, outliers})
\item Scatterplot $\hat{y}$ versus $y$ (\textem{How good is the fit of
the linear model?  Linearity, constant variance}.)
\item Boxplot of residuals ($e_i$) (\textem{normality, outliers})
\end{enumerate}


\shinyGotoTab{name=Examples; text=Goto Examples tab}

\end{shinyTab}

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Examples tab %%%%%%%%%%%%%%%%%%%%%%%%

\begin{shinyTab}{Examples}

\h3{Examples}

\shinyInternalLoc{name=chooseExample}
\shinyDropdown{inputId=exDrop ; label=Choose an example ; 
   choices=blood level, yacht racing, Kittiwake colonies }

%-------- blood level example ----------
\begin{shinyConditionalPanel}{inputId=exDrop; value=blood level}

\h4{Blood level example}

\textbf{Background:}

The data set, called \textem{BloodLevel.dat}, is concerned with the
relationship between the levels of a hormone in the blood (\textsc{
bl.level}) and age(\textsc{age}, measured in weeks) in the first four
weeks of life.  The observations are taken from 25 healthy babies, 5
newborns, 5 in week 1, 5 in week 2, 5 in week 3, and 5 in week 4.
Doctors are interested in describing the relationship between \textsc{
bl.level} and \textsc{age} using a regression model in order to identify
infants whose blood levels of the hormone are unusual compared to
healthy babies.

\begin{shinyCode}
bloodLevel = read.table("BloodLevel.dat", header=TRUE)
bloodModel = lm(Bl.Level ~ Age, data=bloodLevel)
source("http://www.stat.cmu.edu/~hseltman/files/qqn.R")
\end{shinyCode}

\textbf{Data:}

Here are five random lines from the data set:
\begin{shinyCodeOutput}
print(bloodLevel[sample(1:nrow(bloodLevel),5), 1:2], row.names=FALSE)
\end{shinyCodeOutput}

The dependent variable is the blood level of the hormone (\textsc{Bl.Level}), and the
independent variable is the age (\textsc{Age}).  Both are quantitative variables, so
a simple linear regression is one possible reasonable analysis.

\smallskip
\textbf{EDA:}

Here is the univariate EDA:
\begin{verbatim}
summary(bloodLevel$Age)
table(bloodLevel$Age)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(bloodLevel$Age))
print(table(bloodLevel$Age))
\end{shinyCodeOutput}
(Typically table() is reserved for categorical variable, but with only a few
different levels of a quantitative variable, it can be useful.)

\smallskip
\begin{verbatim}
summary(bloodLevel$Bl.Level)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(bloodLevel$Bl.Level))
\end{shinyCodeOutput}


\smallskip
\begin{verbatim}
par(mfrow=c(1,2))
pie(table(bloodLevel$Age), main="Age")
hist(bloodLevel$Bl.Level, xlab="blood level", main="")
\end{verbatim}
\begin{shinyPlot}
par(mfrow=c(1,2))
pie(table(bloodLevel$Age), main="Age")
hist(bloodLevel$Bl.Level, xlab="blood level", main="")
\end{shinyPlot}
(Atypically and due to the small number of different ages studied,
we show a pie chart instead of a histogram for \textsc{Age}.)

\medskip
Here is the bivariate non-graphical EDA (see the above scatterplot for the bivariate graphical EDA):
\begin{verbatim}
with(bloodLevel, cor(Age, Bl.Level))
\end{verbatim}
\begin{shinyCodeOutput}
print(with(bloodLevel, cor(Age, Bl.Level)))
\end{shinyCodeOutput}

%\smallskip
%\begin{verbatim}
%with(bloodLevel, plot(Bl.Level ~ Age))
%\end{verbatim}
%\begin{shinyPlot}
%with(bloodLevel, plot(Bl.Level ~ Age))
%\end{shinyPlot}

\smallskip
At this point we should have a good grasp of the nature of the data, and do
not see evidence of any gross errors in the data.  We can assess the general
relationship between the \textsc{Age} and \textsc{Bl.Level} variables, e.g.,
"older babies tend to have lower blood levels of the hormone".

We are now ready to consider a statistical model.  One model for a linear
relationship between a quantitative IV and a quantitative DV is
"simple linear regression".  For this example, the model would assume that there 
is a linear relationship between the \textsc{Age} and the mean \textsc{Bl.Level},
which seems reasonable here.  It would also assume that at each fixed age, 
the distribution of the blood level of the hormone would follow a normal (Gaussian)
distribution with the same variance at each age.  In the scatterplot, we
see evidence of a different variance at each age, with the highest variances
seen at the youngest ages.  It is nearly impossible for the human eye to
judge the validity of the normality assumption, so we will save that for later
(residual analysis).  

The simple linear regression assumption of independent errors cannot be checked
with the EDA, but we can think about whether we know anything that might put
doubt in that assumption.  If the \textbfem{same} babies were tested at each
age, then the errors (unexplained individual deviations from the group mean at each age)
would probably be similar for the same baby at different ages, violating the assumption.
Or, if the babies interacted with each other in groups, and the nature of the interaction
affected the hormone levels, that might violate the assumption.

For now we will proceed with the regression analysis, with some concern reserved
for possible violation of the equal variance assumption.

\smallskip
\textbf{Regression analysis}
\begin{verbatim}
summary(lm(Bl.Level ~ Age, data=bloodLevel))
\end{verbatim}
\begin{shinyCodeOutput}
  txt = capture.output(summary(bloodModel))[9:18]
  # Fix weird backward quotes that show wrong in Shiny
  signif = grep("Signif", txt)
  badPunct = c("‘|’")
  txt[signif] = gsub(badPunct, "'", txt[signif])
  cat(txt, sep="\n")
\end{shinyCodeOutput}

\smallskip
\textbf{Interpretation:}

The \textbf{intercept estimate}, $b_0=13.3$, is an estimate of the population intercept, $\beta_0$,
which represents the mean of the DV ($Y$, \textsc{Bl.Level}) when the IV equals 0.  In this
case it is the mean hormone blood level for newborns (\textsc{Age}=0).  With a p-value
of $p<<0.0001$, we can reject the null hypothesis that $\beta_0=0$.  But this is
uninteresting, because the null hypothesis is scientifically untenable, and because if the mean
were zero, then the model Gaussian distribution would imply that half of the values
are negative, which is impossible.  We can use the standard error of the coefficient (SE=0.76)
to construct a 95\% confidence interval for $\beta_0$ as coefficient $\pm$ 2SE = [11.7, 14.8].
Then we can conclude that we are 95\% confident that the true mean \textsc{Bl.Level} for newborns
is between those two values.

The \textbf{slope estimate}, $b_1=-2.11$, is an estimate of the population slope, $\beta_1$,
which represents the change in the mean \textsc{Bl.Level} for each 1 unit (year) increase in
\textsc{Age} (starting at any age), the mean of \textsc{Bl.Level} is predicted to fall
by 2.11.  To be specific, we are not claiming that each individual's \textsc{Bl.Level} falls
by 2.11 for each year they age.  Instead we are claiming that if we compute the mean
\textsc{Bl.Level} for any specific age (within the range studied), then we expect that the
mean \textsc{Bl.Level} for all subjects the age one year larger will be 2.11 lower.
This is only a so-called "point estimate" of the mean difference, and does not
include any information on the uncertainty in our estimate caused by estimating it
from a single random sample.  But we can compute an approximate 95\% confidence interval
on the change in mean \textsc{Bl.Level} for each one year rise in age as
$-2.11 \pm 2*0.31 = [-2.73, -1.49]$.  We are 95\% confident that the true yearly
change in mean \textsc{Bl.Level} is between 1.49 and 2.73 hormone units.

We can also interpret the "residual standard error" of 2.20 as an estimate of
$\sigma$, the error standard deviation.  Our best point estimate of the 
s.d. characterizing the spread of \textsc{Bl.Level} values around the mean
for any given age is 2.20.

\smallskip
\h4{Model Checking}
The two standard plots to check model assumptions are the \textbf{residual vs. fit} plot
and the \textbf{quantile normal plot} of residuals.

\begin{shinyPlot}
  plot(residuals(bloodModel)~fitted(bloodModel), xlab="Fitted Values", ylab="Residuals",
       main="Residual vs. Fit for Regression of Bl.Level on Age")
  abline(h=0)
\end{shinyPlot}

\smallskip
The residual vs. fit plot shows the fitted values ($\hat{y_i}$) on the x-axis and
the residual ($e_i$) values on the y-axis.  Points above the y=0 reference line
represent subjects whose observed DV value is higher than predicted based on their IV value,
and vice versa.  If the equal variance assumption is violated, then in some
regions of the plot the vertical spread is larger than in other regions (horizontally).
Here we see greater spread when the fitted (predicted) values are high than when they are low.
If the linearity assumptions is violated then we see either a "smile" or a "frown" pattern
where the vertical mean of the points is not fixed near zero as we move left to right
across the plot.  Here we see a bit of a "smile" pattern.

\smallskip
The quantile normal plot of the residuals is an aid to testing the normality assumption.

\begin{shinyPlot}
  qqn(residuals(bloodModel),
       main="Residual QN Plot for Regression of Bl.Level on Age")
\end{shinyPlot}
Here we see a slight pattern like a side-ways "S" shape, suggesting very mild non-normality
in the form of heavy tails (low residual at the low end, and high residuals at the high end).

\h4{Using transformation to improve the model assumptions:}

When the equal variance and/or normality assumptions are not met, transformation of 
the y variable (DV) to put it on a different scale may fix the problem.  When the
linearity assumption is violated, transformation of the IV, DV, or both may
fix the problem.

We can try to use the log of the \textsc{Bl.Level} which we call \textsc{Log.Bl.Level}.  Note
that it does not matter if we use base 10 logs or natural logs, as long as we keep track
of which one we are using.

Here is the new EDA plot:
\begin{shinyCode}
  bloodLevel[["Log.Bl.Level"]] = log(bloodLevel[["Bl.Level"]])
  bloodModel2 = lm(Log.Bl.Level ~ Age, data=bloodLevel)
\end{shinyCode}
\begin{shinyPlot}
with(bloodLevel, plot(Log.Bl.Level ~ Age))
\end{shinyPlot}
\smallskip

We see that the difference in (vertical) spread from the left to the right
side of the plot is less severe.

\smallskip
Here is the new regression:
\begin{shinyCodeOutput}
  txt = capture.output(summary(bloodModel2))[9:18]
  signif = grep("Signif", txt)
  badPunct = c("‘|’")
  txt[signif] = gsub(badPunct, "'", txt[signif])
  cat(txt, sep="\n")
\end{shinyCodeOutput}

\smallskip
And here is are the residual plots:
\begin{shinyPlot}
  plot(residuals(bloodModel2)~fitted(bloodModel2), xlab="Fitted Values", ylab="Residuals",
       main="Residual vs. Fit for Regression of Log.Bl.Level on Age")
  abline(h=0)
\end{shinyPlot}
\smallskip
We see less evidence of non-linearity and unequal variance.
\begin{shinyPlot}
  qqn(residuals(bloodModel2),
       main="Residual QN Plot for Regression of Bl.Level on Age")
\end{shinyPlot}
\smallskip
Unfortunately, the "heavy tailed" non-normality is a little worse here than in than
original plot.  

\smallskip
\textbf{Conclusion:}  The analysis with \textsc{Bl.Level} on the log scale gives
more reliable conclusions because the violations of the linearity and equal
variance assumptions are smaller.  (Although the normality violation is a bit
larger, the simple linear regression model is quite robust to violations of the
normality assumption.)

Our best fit line is now:
$ \mathrm{log(Bl.Level)} = 2.58 - 0.23 \mathrm{Age}$.  We can express this as
the predicted blood hormone level gets multiplied by $\exp(-0.23)=0.79$ which
is equivalent to a 21\% drop for each additional year of age.

\end{shinyConditionalPanel}


%-------- second example: yacht racing ----------
\begin{shinyConditionalPanel}{inputId=exDrop; value=yacht racing}
\h4{Yacht Racing}

Note:  These data are from the Sydney Morning Herald of 24 December 1994 and
are available at \shinyLinkNew{OZDASL}{http://www.statsci.org/data/oz/sydhob.html}.

\textbf{Background:}

The Sydney-Hobart yacht race starts from Sydney Harbour on Boxing day (December 26) and finishes several days later in Hobart. It is a 630 nautical mile ocean race. The data give the winning times from 1945 to 1993, as they appeared in the Sydney Morning Herald on 24 December, 1994, plus the winning times for 1994 to 1997.

\textbf{Data:}
\begin{shinyCode}
yacht = read.csv("http://www.statsci.org/data/oz/sydhob.txt", sep="\t",
                 as.is=TRUE)[, c("Year", "Time")]
yachtModel = lm(Time ~ Year, data=yacht)
\end{shinyCode}

Here are five random lines from the data set:
\begin{shinyCodeOutput}
print(yacht[sample(1:nrow(yacht), 5), ], row.names=FALSE)
\end{shinyCodeOutput}

The dependent variable is the time to complete the race in minutes (\textsc{Time}), and the
independent variable is year the race was run (\textsc{Year}).  Both are quantitative variables,
so a simple linear regression is one possible reasonable analysis.

\smallskip
\textbf{EDA:}

Here is the univariate EDA:
\begin{verbatim}
summary(yacht$Year)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(yacht$Year))
\end{shinyCodeOutput}


\smallskip
\begin{verbatim}
summary(yacht$Time)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(yacht$Time))
\end{shinyCodeOutput}


\smallskip
\begin{verbatim}
par(mfrow=c(1,2))
hist(yacht$Year, xlab="Year", main="")
hist(yacht$Time, xlab="Time (minutes)", main="")
\end{verbatim}
\begin{shinyPlot}
par(mfrow=c(1,2))
hist(yacht$Year, xlab="Year", main="")
hist(yacht$Time, xlab="Time (minutes)", main="")
\end{shinyPlot}

\medskip
Here is the bivariate non-graphical EDA (see the above scatterplot for the bivariate graphical EDA):
\begin{verbatim}
with(yacht, cor(Year, Time))
\end{verbatim}
\begin{shinyCodeOutput}
print(with(yacht, cor(Year, Time)))
\end{shinyCodeOutput}


\medskip
Here is the bivariate graphical EDA (see the above scatterplot for the bivariate graphical EDA):
\begin{verbatim}
plot(Time~Year, data=yacht, ylab="Time (minutes)")
abline(lm(Time~Year, data=yacht), col=2, lwd=2)
\end{verbatim}
\begin{shinyPlot}
plot(Time~Year, data=yacht, ylab="Time (minutes)")
abline(lm(Time~Year, data=yacht), col=2, lwd=2)
\end{shinyPlot}

%\smallskip
%\begin{verbatim}
%with(yacht, plot(Time ~ Year))
%\end{verbatim}
%\begin{shinyPlot}
%with(yacht, plot(Time ~ Year))
%\end{shinyPlot}

\smallskip
At this point we should have a good grasp of the nature of the data, and do
not see evidence of any gross errors in the data.  We can assess the general
relationship between the \textsc{Time} and \textsc{Year} variables, e.g.,
"races are run faster in later years".

We are now ready to consider a statistical model.  One model for a linear
relationship between a quantitative IV and a quantitative DV is
"simple linear regression".  For this example, the model would assume that there 
is a linear relationship between the \textsc{Time} and \textsc{Year},
which seems reasonable here.  It would also assume that the distribution of the
"errors" (deviations of actual race times from the value predicted by the
regression line) follow a normal (Gaussian)
distribution with the same variance at each year.  In the scatterplot, we
see no evidence of a different variance at different years.
It is nearly impossible for the human eye to
judge the validity of the normality assumption, so we will save that for later
(residual analysis).  

The simple linear regression assumption of independent errors cannot be checked
with the EDA, but we can think about whether we know anything that might put
doubt in that assumption.  If a particularly fast boat won for several years
in a row, then the errors
might be similar for that boat in several years, violating the assumption.

With no initial evidence of assumption violation, we can proceed to the formal linear regression analysis.

\smallskip
\textbf{Regression analysis}
\begin{verbatim}
summary(lm(Time ~ Year, data=yacht))
\end{verbatim}
\begin{shinyCodeOutput}
  txt = capture.output(summary(yachtModel))[9:18]
  # Fix weird backward quotes that show wrong in Shiny
  signif = grep("Signif", txt)
  badPunct = c("‘|’")
  txt[signif] = gsub(badPunct, "'", txt[signif])
  cat(txt, sep="\n")
\end{shinyCodeOutput}

\smallskip
\textbf{Interpretation:}

The \textbf{intercept estimate}, $b_0=144708$ minutes, is an estimate of the population intercept, $\beta_0$,
which represents the mean of the DV ($Y$, \textsc{Time}) extrapolated to where the IV 
(\textsc{Year}) equals 0.  Because year zero is a ridiculous extrapolation, there is no
interpretation of the estimate, and there is no interpretation of the p-value.  But
remember that you do need to use that intercept to make predictions of race time
for any particular year, though predictions too far outside of the observed range of 1945 to 1997
is dangerous (i.e, it could easily lead to incorrect estimates due to lack
of linearity outside of the observed range).

The \textbf{slope estimate}, $b_1=-70.6$, is an estimate of the population slope, $\beta_1$,
which represents the change in the mean \textsc{Time} for each 1 unit (year) increase in
\textsc{Year} (starting at any year), the mean of \textsc{Time} is predicted to fall
by 70.6.

This is only a so-called "point estimate" of the mean difference, and does not
include any information on the uncertainty in our estimate caused by estimating it
from a single random sample.  But we can compute an approximate 95\% confidence interval
on the change in mean \textsc{Time} for each subsequent year as
$-70.6 \pm 2*7.82 = [-82.2, -55.0]$.  We are 95\% confident that the true 
drop in mean \textsc{Time} for each subsequent year 
is between 55.0 and 82.2 minutes.

We can also interpret the "residual standard error" of 870 as an estimate of
$\sigma$, the error standard deviation.  Our best point estimate of the 
s.d. characterizing the spread of \textsc{Times} values around the mean
for all years is 870 minutes.

\smallskip
\h4{Model Checking}
The two standard plots to check model assumptions are the \textbf{residual vs. fit} plot
and the \textbf{quantile normal plot} of residuals.

\begin{shinyPlot}
  plot(residuals(yachtModel)~fitted(yachtModel), xlab="Fitted Values", ylab="Residuals",
       main="Residual vs. Fit for Regression of Minutes on Time")
  abline(h=0)
\end{shinyPlot}

\smallskip
The residual vs. fit plot shows the fitted values ($\hat{y_i}$) on the x-axis and
the residual ($e_i$) values on the y-axis.  Points above the y=0 reference line
represent subjects whose observed DV value is higher than predicted based on their IV value,
and vice versa.  If the equal variance assumption is violated, then in some
regions of the plot the vertical spread is larger than in other regions (horizontally).
Here we see no evidence of a violation of this assumption.
If the linearity assumptions is violated then we see either a "smile" or a "frown" pattern
where the vertical mean of the points is not fixed near zero as we move left to right
across the plot.  Here we see no evidence of a violation of this assumption.

\smallskip
The quantile normal plot of the residuals is an aid to testing the normality assumption.

\begin{shinyPlot}
  qqn(residuals(yachtModel),
       main="Residual QN Plot for Regression of Minutes on Time")
\end{shinyPlot}
Here we see an extremely minor deviation from the reference line in the form  pattern like a side-ways "S" shape, suggesting extremely mild non-normality or, more likely, just due to the random
variation from sample to sample.  Because mild violations of normality have essentially
no effect on regression results, we can proceed as if the assumption were met.


\smallskip
\textbf{Conclusion:}  If the untestable independent errors assumption is not violate,
then this analysis gives
reliable conclusions because there are no detectable violations of the linearity and equal
variance assumptions.

Our best fit line is now:
$ Minutes = 144708 - 70.5 \mathrm{Year}$.  We can express this as
the predicted time will go down by 70.5 minutes for each subsequent year.

Note that it might be a good idea to create a new variable that is "years since
1945" and use this in place of $\mathrm{Year}$ because that would allow
a good meaning for the intercept, namely the predicted race time for 1945.

\end{shinyConditionalPanel}



%-------- third example: Kittiwake Colonies ----------
\begin{shinyConditionalPanel}{inputId=exDrop; value=Kittiwake colonies}
\h4{Kittiwake Colonies}

Note:  These data are from \shinyLinkNew{OzDASL}{http://www.statsci.org/data/general/kittiwak.html}.

\textbf{Background:}

Cairns (\textem{The regulation of seabird colony size: a hinterland model}. \textbf{The American Naturalist}, 134(1), 141-146, 1988) analysed the relation between population and foraging area for seabird colonies. The following table presents their data for 22 black-legged kittiwake (a northern gull) colonies of Scotland's Shetland and Orkney Islands. Area is km2 and Population is the number of breeding pairs. 

\begin{verbatim}
birds = read.table("http://www.statsci.org/data/general/kittiwak.txt", header=TRUE, sep="\t")
\end{verbatim}
\begin{shinyCode}
birds = read.table("http://www.statsci.org/data/general/kittiwak.txt", header=TRUE, sep="\t")
birdsModel = lm(Population ~ Area, data=birds)
\end{shinyCode}
\begin{shinyPlot}
plot(Population ~ Area, data=birds)
abline(lm(Population ~ Area, data=birds), col=2, lwd=2)
\end{shinyPlot}


\textbf{Data:}

Here are five random lines from the data set:
\begin{shinyCodeOutput}
print(birds[sample(1:nrow(birds),5), ], row.names=FALSE)
\end{shinyCodeOutput}

The dependent variable is the population of the colony \textsc{Population}, and the
independent variable is the \textsc{Area} (in $km^2$).  Both are quantitative variables,
so a simple linear regression is one possible reasonable analysis.

\smallskip
\textbf{EDA:}

Here is the univariate EDA:
\begin{verbatim}
summary(birds$Area)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(birds$Area))
\end{shinyCodeOutput}


\smallskip
\begin{verbatim}
summary(birds$Population)
\end{verbatim}
\begin{shinyCodeOutput}
print(summary(birds$Population))
\end{shinyCodeOutput}


\smallskip
\begin{verbatim}
par(mfrow=c(1,2))
hist(birds$Area, xlab="Area (sq. km)", main="")
hist(birds$Population, breaks=10, xlab="Population", main="")
\end{verbatim}
\begin{shinyPlot}
par(mfrow=c(1,2))
hist(birds$Area, xlab="Area (sq. km)", main="")
hist(birds$Population, breaks=10, xlab="Population", main="")
\end{shinyPlot}

\medskip
Here is the bivariate non-graphical EDA (see the above scatterplot for the bivariate graphical EDA):
\begin{verbatim}
with(birds, cor(Area, Population))
\end{verbatim}
\begin{shinyCodeOutput}
print(with(birds, cor(Area, Population)))
\end{shinyCodeOutput}

%\smallskip
%\begin{verbatim}
%with(birds, plot(Population ~ Area))
%\end{verbatim}
%\begin{shinyPlot}
%with(birds, plot(Population ~ Area))
%\end{shinyPlot}

\smallskip
At this point we should have a good grasp of the nature of the data, and do
not see evidence of any gross errors in the data, although there are
a small number of "outliers" with very large area and population.  We can assess the general
relationship between the \textsc{Area} and \textsc{Population} variables, e.g.,
"a larger area corresponds to a higher population".

We are now ready to consider a statistical model.  One model for a linear
relationship between a quantitative IV and a quantitative DV is
"simple linear regression".  For this example, the model would assume that there 
is a linear relationship between the \textsc{area} and the mean \textsc{population}
for islands with that area,
which seems reasonable here.  It would also assume that at each fixed area, 
the distribution of the population would follow a normal (Gaussian)
distribution with the same variance for each area.  In the scatterplot, we
see no evidence of a different variance at each area.
It is nearly impossible for the human eye to
judge the validity of the normality assumption, so we will save that for later
(residual analysis).  

The simple linear regression assumption of independent errors cannot be checked
with the EDA, but we can think about whether we know anything that might put
doubt in that assumption.  If there were, e.g., groups of islands closes enough
for birds to fly between them, this might induce correlation in the errors.

With no initial evidence of assumption violation, we can proceed to the formal linear regression analysis.

\smallskip
\textbf{Regression analysis}
\begin{verbatim}
summary(lm(Population ~ Area, data=birds))
\end{verbatim}
\begin{shinyCodeOutput}
  txt = capture.output(summary(birdsModel))[9:18]
  # Fix weird backward quotes that show wrong in Shiny
  signif = grep("Signif", txt)
  badPunct = c("‘|’")
  txt[signif] = gsub(badPunct, "'", txt[signif])
  cat(txt, sep="\n")
\end{shinyCodeOutput}

\smallskip
\textbf{Interpretation:}

The \textbf{intercept estimate}, $b_0=-735$, is an estimate of the population intercept, $\beta_0$,
which represents the mean of the DV ($Y$, bird \textsc{population}) extrapolated to where the IV 
(\textsc{area}) equals 0.  This is a meaningless extrapolation.
(Note that possible word confusion in this example.  "Population intercept"
refers to a theoretical population of all possible bird colonies as opposed to the
sample intercept which refers to the colonies in this study.  In contast, "bird population"
refers to the number of birds in a particular colony.)

The \textbf{slope estimate}, $b_1=3.30$, is an estimate of the population slope, $\beta_1$,
which represents the change in the mean bird \textsc{population} for each 1 unit increase in
\textsc{area} (starting at any area), the mean of the bird \textsc{population} is predicted to rise by 3.30.  
This is only a so-called "point estimate" of the mean difference, and does not
include any information on the uncertainty in our estimate caused by estimating it
from a single random sample.  But we can compute an approximate 95\% confidence interval
on the change in mean bird \textsc{population} for each one unit increase in area as
$3.30 \pm 2*0.58 = [2.14, 4.46]$.  We are 95\% confident that the true 
drop in mean bird \textsc{populatoin} for each one unit increase in the are 
is between 2.14 and 4.46.

We can also interpret the "residual standard error" of 2579 as an estimate of
$\sigma$, the error standard deviation.  Our best point estimate of the 
s.d. characterizing the spread of bird \textsc{population} values around the mean
for any given area is 2579.

\smallskip
\h4{Model Checking}
The two standard plots to check model assumptions are the \textbf{residual vs. fit} plot
and the \textbf{quantile normal plot} of residuals.

\begin{shinyPlot}
  plot(residuals(birdsModel)~fitted(birdsModel), xlab="Fitted Values", ylab="Residuals",
       main="Residual vs. Fit for Regression of response on amount")
  abline(h=0)
\end{shinyPlot}

\smallskip
The residual vs. fit plot shows the fitted values ($\hat{y_i}$) on the x-axis and
the residual ($e_i$) values on the y-axis.  Points above the y=0 reference line
represent subjects whose observed DV value is higher than predicted based on their IV value,
and vice versa.  If the equal variance assumption is violated, then in some
regions of the plot the vertical spread is larger than in other regions (horizontally).
We do see evidence of unequal variance because the vertical spread of the data points
is smaller on the left side of the plot (low fitted values, which in this case, with 
a single IV and a positive $\hat{\beta_1}$, corresponds to lower amounts of the
substance),  than on the right.
If the linearity assumptions is violated then we see either a "smile" or a "frown" pattern
where the vertical mean of the points is not fixed near zero as we move left to right
across the plot.  Here we see a "smile" when we connect the vertical center of each band
from left to right, which is moderately strong evidence of violation of the linearity assumption.

This means that the above interpretations of intercept and slope are inappropriate.


\smallskip
The quantile normal plot of the residuals is an aid to testing the normality assumption.
It is unreliable in the case of strong non-linearity, so we skip that check for now.


\smallskip
\h4{Using transformation to improve the model assumptions:}

When the equal variance and/or normality assumptions are not met, transformation of 
the y variable (DV) to put it on a different scale may fix the problem.  When the
linearity assumption is violated, transformation of the IV, DV, or both may
fix the problem.

One possibility is to use the log of the \textsc{population} which we call
\textsc{logPopulation}.  Note
that it does not matter if we use base 10 logs or natural logs, as long as we keep track
of which one we are using.

Here is the new EDA plot:
\begin{shinyCode}
  birds$logPopulation = log(birds$Population)
  birdsModel2 = lm(logPopulation ~ Area, data=birds)
\end{shinyCode}
\begin{shinyPlot}
with(birds, plot(logPopulation ~ Area))
birdsModel2 = lm(logPopulation ~ Area, data=birds)
abline(birdsModel2, col=2, lwd=2)
\end{shinyPlot}

\smallskip
Here is the residual vs. fitted plot:

\begin{shinyPlot}
modelbirds2 = lm(logPopulation ~ Area, birds)
plot(residuals(modelbirds2) ~ fitted(modelbirds2), xlab="fitted", ylab="residual",
     main="Simple regression of log(response) on log(amount)")
abline(h=0)
\end{shinyPlot}

This looks fine, with a random scatter around the horizontal line at zero and 
the same vertical spread at the left, middle and right portions of the plot.

\smallskip
Now we check the quantile-normal plot:
\begin{shinyPlot}
  qqn(residuals(birdsModel2),
       main="Residual QN Plot for Regression of log(Population) on Area")
\end{shinyPlot}

We conclude that the residuals are reasonable close to Normal in distribution.

\smallskip
Now we can re-fit the model using logPopulation.
\begin{verbatim}
summary(lm(logPopulation ~ Area, data=birds))
\end{verbatim}
\begin{shinyCodeOutput}
  txt = capture.output(summary(birdsModel2))[9:18]
  # Fix weird backward quotes that show wrong in Shiny
  signif = grep("Signif", txt)
  badPunct = c("‘|’")
  txt[signif] = gsub(badPunct, "'", txt[signif])
  cat(txt, sep="\n")
\end{shinyCodeOutput}

We now have the interpretation that the natural log of the population goes up
on average by 0.000910 for every 1 $km^2$ increase in area.  Because a one
$km^2$ change in area is miniscule, it would be better to look at, say,
a 1000 $km^2$ change to get an interpretable explanation: the natural log of the bird population
goes up on average by 0.910 for every 1000 $km^2$ increase in area.  Then we
note that $e^0.910 = 2.48$, and we can conclude that for each 1000 $km^2$ increase
in area the bird population is \textbf{multiplied} by 2.48 (on average).

\smallskip
\textbf{Conclusion:}  By careful examination of the residual vs. fit plot,
we saw that we needed to modify our model to use the log of the population
as our DV to achieve a linear relationship.  Our new model looks like a good one.


\end{shinyConditionalPanel}


\shinyInternalLink{name=chooseExample; text=Choose another example}

\smallskip
\shinyGotoTab{name=Post-test; text=Goto Post-test tab}

\end{shinyTab}

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Post-test tab %%%%%%%%%%%%%%%%%%%%%%%%

\begin{shinyTab}{Post-test}

\h3{Post-test}

\shinyQuiz{bank=postQuiz; addNumbers=TRUE; randomizeQuestions=FALSE}

\smallskip
\h4{Congratulations: you have completed the simple linear regression review app.}

\end{shinyTab}

\end{document}
