MATH 6630: Applied Statistics I
Fall 2023
Where there is no uncertainty, there cannot be truth – Richard Feynman
Negative of Himmelblau’s function:
\(f(x,y) = - (x^2 + y -11)^2 - (x + y^2
-7)^2\), see p. 57, problem 2.7 of the textbook.
Announcements
September 18:
- R tutorial at 3:30 on Wednesdays starting September 20. On Zoom at https://yorku.zoom.us/j/92413559596
- General tutorial at 3:30 on Fridays at https://yorku.zoom.us/j/92413559596
- ‘Discussion hour’ idea cancelled, the tutorial will just keep going
until no one is left!
- We’ll have an ‘office hour’ determined later.
September 14:
- See sample quiz questions. Use the link below, in ‘Quick links’.
- Remember the first tutorial and ‘discussion hour’, on Friday,
September 15:
- Tutorial at 3:30 at https://yorku.zoom.us/j/92413559596
Discussion hour at 4:30 at https://yorku.zoom.us/my/georgesmonette
Quick links
- Current date
- Course description
- Piazza class forum
- Optional tutorial zoom link, Friday 3:30 until the last person leaves. I might have to wrap up around 5 or 5:30.
- R focused tutorial zoom link, Wednesdays at 3:30 - 4:30. These will continue while the demand lasts.
Discussion hour, Friday 4:30cancelled, just an extension of the tutorial instead.- “Office hour == ‘drop-in’ hour” but with possibility of privacy if necessary, will be announced.
- Eclass – only for grades
- Recordings
- The sound is awful in the September 13 recording but it will be much better in future recordings.
- The recording for the tutorial on September 15 failed, but I’m confident it will work in the future.
- Be warned that in the September 22 tutorial video, I make a mistake at 18:52, which Evan makes me realize at 19:43. So don’t get confused by the discussion between those two time points!
- Quizzes
- Sample quiz questions
- Questions
Calendar
Classes meet on Mondays and Wednesdays from 1:00 to
2:30. The location on Mondays is Ross South 123 and, on Wednesdays, Life
Sciences Building (LSB) 107.
Optional Tutorial: Fridays from 3:30 pm to 4:30 pm on
Zoom at Tutorial Zoom
link
- The tutorial scheduled for Friday, September 29 will take place on Monday, October 2 at 3:30 pm.
Informal chat: Tentatively Fridays, 4:30 pm to 5:30
pm at Office hour Zoom
link
Instructor: Georges Monette 1 2
Important dates: 2023-2024 Academic
Calendar
Textbook: Givens, G.H, Hoeting J.A.
(2013) Computational Statistics 2nd ed., Wiley. Some copies
available in bookstore. Errata for the
second edition
Email: Messages about course content should be posted
publicly to Piazza so the discussion can engage and benefit everyone.
You may post messages and questions as personal messages to the
instructor. If they are of general interest and don’t contain personal
information, they will usually be made public for the benefit of the
entire class unless you specifically request that the message remain
private.
Day 1: Wednesday, September 6
Files used today: Regression Review
Assignment 1 (individual): Setting things up
Due Wednesday, September 13, 9 pm
- Summary:
- Install R and RStudio
- Get a free Github account
- Install git on your computer
- Connect with Piazza, create a LOG post and introduce yourself
- 1. Install R and RStudio following these instructions. If you already have R and RStudio, update them to the latest versions.
- 2. Get a free Github account: If you don’t have
one, follow this
excellent source of advice from Jenny Bryan on using git and Github
with R and RStudio.
- Be sure to read Jenny Bryan’s advice before choosing a name.
- CAUTION: Avoid installing ‘Github for Windows’ (which is not the same as ’Git for Windows) from the Github site at this stage.
- An excellent additional reference is Hadley Wickham’s chapter on Git and Github in his book on R Packages
- 3. Install git on your computer using these instructions.
- If you are curious about how git is used have a look at this tutorial! You will be using it mainly through RStudio and won’t need to master how to use with a command-line interface.
- As a final step: In the RStudio menu click on
Tools | Global Options ... | Terminal. Then click on the box in Shell paragraph to the right of New terminals open with:- On a PC, select Git Bash
- On a Mac, select Bash
- On a PC, select Git Bash
- You don’t need to do anything else at this time. We will see how to set up SSH keys to connect to Github through the RStudio terminal in a few lectures.
- 4. Connect with Piazza and post about your
experiences installing the above:
- Join the MATH 6630 Piazza site by going to this URL: piazza.com/yorku.ca/winter2023/math6630
- Use the access code blackwell when prompted.
- Create a post on Piazza to introduce yourself to your colleagues:
- The title should read:
Introduction:followed by the name you prefer to use in class, e.g.Jon Smith - The first line should show the name of your Github account,
e.g.
Github: jsmith. If you don’t have one already, leave this blank for now. You can come back and edit it later.
- Follow this with information on what computing languages you are familiar with? Which one are you proficient with?
- Share some interests: hobbies, favourite musicians, movies, restaurants, etc.
- Click on the
assn1folders when you submit your post.
- The title should read:
- Create another post entitled ‘Getting started’.
- Describe problems you encountered installing R, RStudio, and git
- Describe problems registering with Github
- How could the instructions be improved to make the process easier?
- If you couldn’t complete the installation, describe the problem or error message(s) you encountered.
- Click on the on the
assn1folder before submitting your post. - If you are hoping to get an answer to your post, select
Questionas thePost Type.
- Join the MATH 6630 Piazza site by going to this URL: piazza.com/yorku.ca/winter2023/math6630
Assignment 2 (teams) \(p\)-values
- Deadlines: See the course description for the meaning of these
deadlines.
- Wednesday, September 20 at noon
- Friday, September 22 at 9 pm
- Sunday, September 24 at 9 pm
This exercise will be done in teams of 4 or 5 members which will be posted on Piazza. All members of a team work on a similar problem but with different parameters. This way, you can help each other with the concepts but each team member must, ultimately, do their own work.
The purpose of the assignment is to explore the meaning of p-values. Before starting, stop and reflect on what it means for an experiment to ‘achieve’ a p-value of 0.049. What meaning can we give to the quantity ‘0.049’? How is it related to the probability that the null hypothesis is correct?
To keep things very simple suppose you want to test \(H_0: \mu =0\) versus \(H_1: \mu > 0\) (note: the problem is also easy but takes a bit more coding if you assume a two-sided alternative: \(H_0: \mu =0\) versus \(H_1: \mu \neq 0\) – to help focus on the essential concepts, we’ll use a one-sided alternative) and you are designing an experiment in which you plan to take a sample of independent random variables, \(X_1, X_2, ... , X_n\) which are iid \(\textrm{N}(\mu,1)\), i.e. the variance is known to be equal to 1. You plan to use the usual test based on \(\bar{X}_n\) rejecting \(H_0\) for values of \(\bar{X}_n\) that are far from 0.
An applied example would be testing for a change in value of a response when all subjects are submitted to the same conditions and the measurement error of the response is known. In that example \(X_i\) would be the ‘gain score’, i.e. post-test response minus the pre-test response exhibited by the \(i\)th subject.
Note: Do the assignment using R and using a Markdown file. You can use the Markdown sample file in Assignment 1 as a template. Upload both a pdf file and your R script to Piazza.
Let the selected probability of Type I error be \(\alpha = 0.05\). Depending on your team, you have a budget that allows you to collect a sample of \(n\) observations as follows:
| Team | \(n\) |
|---|---|
| A: Akaike (Hirotsugu Akaike) | 2 |
| B: Blackwell (David Blackwell) | 3 |
| C: Cox (Gertrude Mary Cox) | 5 |
| D: David (Florence Nightingale David) | 10 |
| E: Efron (Bradley Efron) | 15 |
| F: Freedman (David Freedman) | 20 |
Depending on your position on the team, you are interested in answering the questions below using the following values of \(\mu_1\):
| Position in Team | \(\mu_1\) | Cohen’s terminology for effect size: \(\mu_1/\sigma\) |
|---|---|---|
| 1 | all | combine ROC curves for results below |
| 2 | 0.2 | small effect size |
| 3 | 0.5 | medium effect size |
| 4 | 0.8 | large effect size |
| 5 | 1 | very large effect size |
- What is the probability that \(p \le 0.05\) if \(H_0: \mu = 0\) is true?
- What is the probability that \(p \le 0.05\) if \(\mu = \mu_1\)?
- What is the power of this test if \(\mu = \mu_1\)?
- Plot the ROC curve: Power (\(\beta\) as a function of \(\alpha\)). To generate the ROC curve, use a reasonable range of values for critical values of the test and plot the probability of rejecting under \(H_0\) in the horizontal direction versus the probability of rejecting of under \(H_A\) in the vertical direction.
- What is the interpretation of the area under the ROC curve?
- Suppose that you collect the data and that the observed \(p\)-value is 0.049. What can you say about the probability that \(H_0\) is true?
- Suppose that, before running the experiment, you were willing to
give \(H_0\) and \(H_1: \mu = \mu_1\) equal probability of
1/2.
- What is the probability that \(H_0\) is true given the event that \(p \le 0.05\)?
- What is the probability that \(H_0\) is true given the event that \(p = 0.049\)?
- What is the probability that \(H_0\) is true given the event that \(p = 0.005\)?
- What is the probability that \(H_0\) is true given the event that \(p = 0.0005\)?
- What is the probability that \(H_0\) is true given the event that \(p \le 0.05\)?
- As a team discuss and post your reflections: Hypothesis testing is often presented as a process that parallels that of determining guilt in a criminal process. We start with a presumption of innocence, i.e. that \(H_0\) is true. We then hear evidence and consider whether it contradicts the presumption of innocence ‘beyond a reasonable doubt.’ Suppose we quantify the presumption of innocence to mean that the prior probability of guilt is small, e.g. \(P(H_1) \le .05\). How small an observed \(p\)-value do you need to obtain in order to ‘flip’ the presumption of innocence to ‘guilt beyond a reasonable doubt’ if that is defined as \(P(H_0 | \mathrm{data}) \le .05\).
- Ditto: What \(p\)-value would we need if the presumption of innocence and guilt beyond a reasonable doubt correspond to \(P(H_1) \le 0.001\) and \(P(H_0|\mathrm{data}) \le 0.001\), respectively?
- Ditto: Courts have often adopted a criterion of \(p < 0.05\) in imitation of the common practice among many researchers. Comment on the possible consequences.
- Ditto: Have a look at this xkcd cartoon. How does the Bayesian statistician derive a probability to make a decision in this example? Show the details.
To delve more deeply into these issues you can read Wasserstein and Lazar (2016) and Wasserstein, Schirm, and Lazar (2019). Concerns about \(p\)-values have been around for a long time, see Schervish (1996). For a short overview see Bergland (2019). There are two interesting and influential papers by John Ioannidis (2005), (2019).
For an entertaining take on related issues see Last Week Tonight: Scientific Studies by John Oliver (2016) (warning: contains strong language and political irony that many could consider offensive – watch at your own risk!).
Preparing for the next class: Play your way through these scripts:
Step your way through the scripts and experiment as you go along. If you have questions or comments, post them on Piazza.
These scripts give you a tour of various features of R. They are meant to provide reminders of things you can do when a need arises.
What’s the best way to learn R? That depends on your learning style, your previous experience and how deep you want to go.
See the discussion in ‘Day 2’. Start or join a discussion on Piazza on ‘best ways to learn R’.
Day 2: Monday, September 11
Correction: There was a problem installing the ‘spida2’ package in R. The workaround is to install the ‘gplots’ package before installing the ‘spida2’ package:
install.packages('gplots')
devtools::install_github('gmonette/spida2')
The ‘spida2’ package should be reinstalled regularly since it’s likely to be updated frequently during the course. Once ‘gplots’ has been installed once, it does not have to be reinstalled each time you reinstall ‘spida2’.
Tutorials: Start this week on Friday, September 15 from 3:30 pm to 4:30 pm.
Student hours: Start this week on Friday, September 15, from 4:30 pm to 5:30 pm.
Teams: I will assign teams this evening for those who have already registered with Piazza so you can start working on assignments. So register soon if you want to have a team and, as importantly, your position on your team.
Learning R
Why R? What about SAS, SPSS, Python, Julia, among others?
SAS is a very high quality, intensely engineered, environment for statistical analysis. It is widely used by large corporations. New procedures in SAS are developed and thoroughly tested by a team of 1,000 or more SAS engineers before being released. It currently has more than 300 procedures.
R is an environment for statistical programming and development that
has accumulated many somewhat inconsistent layers developed over time by
people of varying abilities, many of whom work largely independently of
each other. There is no centralized quality testing except to check
whether code and documentation run before a new package is added to R’s
main repository, CRAN. When this page was last updated (at
format(Sys.time(), '%H:%M') on
format(Sys.time(), '%B %d %Y')), CRAN had 20,119
packages.
In addition, a large number of packages under development are available through other repositories, such as Github.
The development of SAS began in 1966 and that of R (in the form of its predecessor, S, at Bell Labs) in 1976.
The ‘design’ of R (using ‘R’ to refer to both R and to S) owes a lot to the design of Unix. The idea is to create a toolbox of simple tools that you link together yourself to perform an analysis. Unix, now mainly as Linux, 3 commands were simple tools linked together by pipes so the output of one command is the input of the next. To do anything you need to put the commands together yourself.
The same is true of R. It’s extremely flexible but at the cost of requiring you to know what you want to do and to be able to use its many tools in combination with each other to achieve your goal. Many decisions in R’s design were intended to make it easy to use interactively. Often the result is a language that is very quirky for programming.
SAS, in contrast, requires you to select options to run large procedures that purport to do the entire job.
This is an old joke: If someone publishes a journal article about a new statistical method, it might be added to SAS in 5 to 10 years. It won’t be added to SPSS until 5 years after there’s a textbook written about it, maybe another 10 to 15 years after its appearance in SAS.
It was added to R two years ago because the new method was developed as a package in R long before being published.
So why become a statistician? So you can have the breadth and depth of understanding that someone needs to apply the latest statistical ideas with the intelligence and discernment to use them effectively.
So expect to have a symbiotic relationship with R. You need R to have access to the tools that implement the latest ideas in statistics. R needs you because it takes people like you to use R effectively.
The role of R in this course is to help us
- have access to tools to expand our ability to explore and analyze data, and
- learn how to develop and implement new statistical methods. i.e. learn how to build new tools
- deepen our understanding of the use of statistics for scientific discovery as well as for business applications
It’s very challenging to find a good way to ‘learn R’. It depend on where you are and where you want to go. Now, there’s a plethora of on-line courses. See the blog post: The 5 Most Effective Ways to Learn R
In my opinion, ultimately, the best way is to
- play your way through the ‘official’ manuals on CRAN starting with ‘An Introduction to R’ along with ‘R Data Import/Export’. Note however that these materials were developed before the current mounting concern with reproducible research and some of the advice should be deprecated, e.g. using ‘attach’ and ‘detach’ with data.frames.
- read the CRAN task views in areas that interest you.
- Have a look at the over 496,438 questions tagged ‘r’ on stackoverflow.
- At every opportunity, use R Markdown documents (like the sample script you ran when you installed R) to work on assignments, project, etc.
Using R is like playing the piano. You can read and learn all the theory you want, ultimately you learn by playing.
We will discuss the scripts CAR_1.R and CAR_2.R in the tutorial this week. Copy them to a file in RStudio and play with them line by line. You can also run an entire script with Control-Shift-K to see how scripts can be used for reproducible research.
Exercises:
- From New 6630
questions
- 5.1, 5.2, 5.3, 5.4, 5.5,
- 5.6.23, 5.6.24, 5.6.25, 5.6.26, 5.6.27,
- 6.1, 6.2, 6.3, 6.4, 6.5,
- 7.4, 7.5, 7.6, 7.7, 7.8,
- 8.1, 8.2, 8.3, 8.4, 8.5,
- 8.6, 8.7, 8.8, 8.9, 8.10,
- 8.18.a, 8.18.b, 8.18.c, 8.18.d, 8.18.e,
- 8.36.a, 8.36.b, 8.36.c, 8.36.d, 8.36.e,
- 8.51.a, 8.51.b, 8.51.c, (write functions that would work on matrices of any size), 8.61.a, 8.62.a
- 12.1, 12.3, 12.5, 12.7, 12.9,
Assignment 3 (teams)
- Do the exercises above. Each member of a team does 10 questions. Use
the random permutation of the numbers 1 to 5 below to determine the
order of the starting exercise for each member of the team. The first
member does the question whose order is the first random number, the
second student does the question whose order is the second random number
and so on. Then you continue to do the all of
the fifth succeeding questions. For example, the third student in each team gets the random number 2 in the sequence below and would do the 2nd, 7th, 12th, 17th, 22nd, 27th, 32nd, 37th, 42nd and 47th questions. - In other words: Trying to make things less
complicated! Note that each of the 10 lines of questions above has 5
questions. So
- Team Member #1 uses the first random number and does the 3rd question in each line,
- Team Member #2 does the 1st question in each line,
- Team Member #3 does the 2nd question in each line,
- Team Member #4 does the 4th question in each line, and
- Team Member #5 does the 5th question in each line.
- Random numbers: 3 1 2 4 5
- Deadlines: See the course description for the meaning of these
deadlines.
- Friday, September 29 at noon
- Sunday, October 1 at noon
- Monday, October 2 at 9 pm
- IMPORTANT:
- Upload the answer to each question in a single Piazza post (post it
as a Piazza ‘Note’, not as a Piazza ‘Question’) with the title: “A3 5.1”
for the first question, etc.
- You can answer the question directly in the posting or by uploading a pdf file and the R script or Rmarkdown script that generated it.
- When providing help or comments, do so as “followup discussion”.
- Upload the answer to each question in a single Piazza post (post it
as a Piazza ‘Note’, not as a Piazza ‘Question’) with the title: “A3 5.1”
for the first question, etc.
If you want a deeper understanding of R, consider reading Hadley Wickham, Advanced R
Day 3: Wednesday, September 13
Coninuation of Day 2
Day 4: Monday, September 18
Files used today:
- Regression Review
- Leverage and Mahalanobis Distance:
- Play through the Simple Regression R script to get information on using Added Variable Plots (also known as Partial Regression Leverage Plots) and Partial Residual Plots (also known as Component-Plus-Residual Plots).
- The script has a short section on ROC curves at the end.
- Continue working on your own with sample scripts in R/CAR_1 and R/CAR_2.
Comments on Leverage and Mahalanobis Distance:
Some formulas are worth remembering because they work for multiple regression as well as simple regression. We’ve seen: \[ SE(\hat{\beta}_i) = \frac{1}{\sqrt{n}}\frac{s_e}{s_{x_i|others}} \] Here’s another one: You know that the diagonal elements, \(h_{ii}\) of the hat matrix \[ H = X(X'X)^{-1}X' \] are the leverages showing how much a change in one particular \(y_i\) would change the fitted value \(\hat{y}_i\). i.e. \[ h_{ii} = \frac{\partial \hat{y}_i}{\partial y_i} \] It turns out that there’s nothing mysterious about where the \(h_{ii}\)s come from. They have a simple relationship with the squared Mahalanobis distance in the predictor space: \[ h_{ii} = \frac{1}{n}(1 + Z_i^2) \] where, in simple regression, \(Z_i\) is the Z-score (standardized with the MLE estimate of standard deviation, dividing by \(n\) instead of \(n-1\)). \(|Z_i|\) also happens to be Mahalanobis distance in one dimension.
The formula generalizes to multiple regression if you take \(Z_i\) to be Mahalanobis distance, so Mahalanobis distance is just a generalization of the Z-score to higher dimensions.
See the links above on Leverage and Mahalanobis Distance.
Comments:
- Start reading Hadley Wickham, Advanced R Chapters 1-3 and work on the exercises.
- Post questions about R on Piazza or come to the tutorials.
Comment on AVPs: The fact that a simple regression applied to the AVP leads to an estimate of the slope that
- is the same as that of the associated multiple regression,
- has a correlation equal to the partial correlation of the predictor in the multiple regression,
- has the same residuals as the residuals in the multiple regression,
- has the same residual standard error as that of the multiple regression, and finally
- has the same standard error for the estimation of the slope as that of the multiple regression
provided one calculates these quantities using maximum likelihood estimates (it’s all asymptotically true if one uses the degree-of-freedom adjusted estimates), is the content of the Frisch-Waugh-Lovell Theorem that was ‘discovered’ in steps over a period of approximately 30 years. The first version applied to lagged time series and was discovered by Ragnar Frisch. He got the first Nobel Prize in Economics in 1931, for his general contributions, not just the start of this theorem. However, like other statistical ‘objects’ that have names associated with them, the actual original theorem was published by Udny Yule in 1911 (I am shamelessly using Wikipedia, and its gifted contributors, as a source). The theorem was refined over time by Frederick Waugh and Michael Lovell who each added portions of it. Lovell is known for an elegant proof.
The theorem shows how the AVP is a window into the way a multiple regression determines an individual coefficient.
Possible quiz questions:
Look at the link to possible questions for the quiz on Wednesday under “Quick Links” above.
Simpsons’ Paradox:
- The principle behind Simpsons’ Paradox is the same whether the variables involved are continuous or categorical. The many numerical examples don’t give me a feeling of understanding what’s behind it. I need a way to visualize it. But visualizations are quite different depending on the character of each of the X, Y and Z variables.
- Perhaps the easiest case to visualize involves a continuous X and Y and a categorical conditioning variable Z.
- See this example I made
up.
- What might be plausible names for the X and Y variables?
- Note: I found a graphic like this somewhere on the web but I couldn’t find it again for this class, so this is my amateurish attempt at recreating it. If you find something similar on the web, please post the link to Piazza.
- With apologies to anyone who has never seen “The Simpsons”.
- The rotating 3d plot at the top of this page is an example in which all three variables are continuous.
- We’ll see some interesting examples where all three variables are dichotomous, using Paik-Agresti diagrams.
Day 5: Wednesday, September 20
Links:
- More references on R:
- Introduction
to R
- very thorough discussion of basic topics in R
- Working with Data
- Regression in R
- Introduction
to R
Day 6: Monday, September 25
Links:
Day 7: Wednesday, September 27
Announcement:
- The tutorial scheduled for Friday, September 29 will take place on Monday, October 2 at 3:30 pm.
Links:
Background material:
- Causal
Questions
- Overview of confounders, mediators, colliders and moderators
- Paik-Agresti diagram
- Lord’s Paradox
- Understanding
Output in Regression
- This is a bit old and needs updating but it might be useful when you explore regression models and want to ask meaningful creative questions.
- Three
Basic Theorems
- These are proofs of things we use in the course. You are not responsible for these proofs but the techniques, primarily using the properties of orthogonal projections (i.e. symmetric idempotent matrices) are very useful when working with the theory of linear regression.
- The second theorem is the proof that the AVP does what we claim it does. It’s the ‘Frisch-Waugh-Lovell’ theorem of econometrics.
- The third theorem shows why linear propensity scores work because
they produce the same \(\hat{\beta}\)
as a multiple regression with the same predictors as those use for the
propensity score. This serves as a warning that the choice of predictors
for the propensity score isn’t different from the choice of predictors
for a multiple regression model: i.e. you still need to block backdoor
paths and avoid mediators and colliders.
- The Causal
Zoo
- Simulation based on a causal graph showing the consequences of using
models that include or exclude various types of variables. See the graph
on page 9 summarizes the consequences of different choices by plotting
the expected value of \(\hat{\beta}_X\)
versus its variance. The ‘causal effect’ of \(X\) is 4, so any model whose estimator has
a different expected value is biased. Look at what happens if you
include the collider in this model:
Y ~ Z + Col. This shows how models matter when you want to give a coefficient a causal interpretation.
- Simulation based on a causal graph showing the consequences of using
models that include or exclude various types of variables. See the graph
on page 9 summarizes the consequences of different choices by plotting
the expected value of \(\hat{\beta}_X\)
versus its variance. The ‘causal effect’ of \(X\) is 4, so any model whose estimator has
a different expected value is biased. Look at what happens if you
include the collider in this model:
Assignment 4 (teams)
- Question
16 entitled ‘Paradoxes, Fallacies and One Correct Statement’ has 23
statements. Your assignment is to discuss and comment on each of these
statements. Following the usual process:
- team member 1 does question 3 and every fifth question after that,
- team member 2 does question 5 and every fifth question after that,
- team member 3 does question 1 and every fifth question after that,
- team member 4 does question 2 and every fifth question after that,
- team member 5 does question 4 and every fifth question after that.
- Deadlines:
- Friday, October 6
- Sunday, October 15 (there’s lots of time to discuss these questions)
- Monday, October 16
Day 8: Monday, October 2
Announcement:
- The tutorial scheduled for Friday, September 29 will take place today at 3:30 pm at https://yorku.zoom.us/j/92413559596
Sample questions for the quiz on Wednesday:
- “If you drop a variable that is not significant in a multiple linear
regression, the other estimated coefficients should not change very
much.”
Comment.
If this statement is false describe how you would look for a counterexample. If this statement is false provide a plausible counterexample. - “If a linear regression on two variables and an intercept leads to
the rejection of the joint null hypothesis: \(H_0: \beta_1 = \beta_2 = 0\), then it must
be possible to reject at least one of the hypotheses: \(H_0: \beta_1 = 0\) or \(H_0: \beta_2 = 0\).”
Comment. If this statement is false describe how you would look for a counterexample. If this statement is false provide a plausible counterexample. - Consider a linear regression on two variables and an intercept with
coefficients \(\beta_1\) and \(\beta_2\). The simple linear regressions
have regression coefficients \(\gamma_1\) and \(\gamma_2\). Is it possible for a data set
to lead to rejection of the hypothesis \(H_0:
\beta_1 = \beta_2 = 0\) yet not reject either \(H_0: \gamma_1 = 0\) nor \(H_0: \gamma_2 = 0\).
If so, provide a plausible example. What does this say about the ability of forward stepwise regression to find a model that fits well. - In multiple regression, when would you want to exclude a variable even though it is highly significant? Give an example.
- In multiple regression, when would you want to include a variable even though it is not significant? Give an example.
Links:
Background material:
- Causal
Questions
- Overview of confounders, mediators, colliders and moderators
- Paik-Agresti diagram
- Lord’s Paradox
- Understanding
Output in Regression
- This is a bit old and needs updating but it might be useful when you explore regression models and want to ask meaningful creative questions.
- Three
Basic Theorems – slightly revised
- See also Three Basic Theorems Notes
- These are proofs of things we use in the course. You are not responsible for these proofs but the techniques, primarily using the properties of orthogonal projections (i.e. symmetric idempotent matrices) are very useful when working with the theory of linear regression.
- The second theorem is the proof that the AVP does what we claim it does. It’s the ‘Frisch-Waugh-Lovell’ theorem of econometrics.
- The third theorem shows why linear propensity scores work because
they produce the same \(\hat{\beta}\)
as a multiple regression with the same predictors as those use for the
propensity score. This serves as a warning that the choice of predictors
for the propensity score isn’t different from the choice of predictors
for a multiple regression model: i.e. you still need to block backdoor
paths and avoid mediators and colliders.
- The Causal
Zoo
- Simulation based on a causal graph showing the consequences of using
models that include or exclude various types of variables. See the graph
on page 9 summarizes the consequences of different choices by plotting
the expected value of \(\hat{\beta}_X\)
versus its variance. The ‘causal effect’ of \(X\) is 4, so any model whose estimator has
a different expected value is biased. Look at what happens if you
include the collider in this model:
Y ~ Z + Col. This shows how models matter when you want to give a coefficient a causal interpretation.
- Simulation based on a causal graph showing the consequences of using
models that include or exclude various types of variables. See the graph
on page 9 summarizes the consequences of different choices by plotting
the expected value of \(\hat{\beta}_X\)
versus its variance. The ‘causal effect’ of \(X\) is 4, so any model whose estimator has
a different expected value is biased. Look at what happens if you
include the collider in this model:
Day 9: Wednesday, October 4
Continuation of previous
Day 10: Monday, October 16
Review
Day 11: Wednesday, October 18
Midterm test
Day 12: Monday, October 23
Links:
- Computational Statistics data sets and scripts
- Slides for Chapter 1
- Chapter 1 Exercises:
- Work out the details of the score moments for a sample of size \(n\) from the \(N(\mu,\sigma)\) distribution.
- Find the canonical parameters, canonical sufficient statistics and cumulant-generating function for a sample from the \(N(\mu, \sigma)\) distribution treated as an exponential family.
- Derive the posterior density for \(\lambda\) given an observation from an exponential(\(\lambda\)) distribution using an exponential(\(\theta\)) prior for \(\lambda\). Compare with a) a posterior density using an improper uniform prior on \(\lambda\). Compare the modes, expected values and the standard deviations of the posteriors.
- Find the profile likelihood for \(\mu\) and for \(\sigma\) for a sample of \(n\) from the \(N(\mu,\sigma)\).
- Find the estimated variance matrix for a Wald test based on the Hessian evaluated at the MLE for a sample of \(n\) from the \(N(\mu,\sigma)\) distribution.
Day 13: Wednesday, October 25
Links:
Assignment 5 (teams)
- See the Piazza Team post to see which of the exercises listed under the previous day and which problems in Chapter 2 of the textbook should be done by each team member.
- Deadlines (updated):
- Friday, November 10
- Sunday, November 12
- Monday, November 13
Day 13: Monday, October 30
Links:
Some sample questions for quiz on November 1:
Note: You should know the definitions of the Bernoulli, Binominal, Multinomial, Normal, Multivariate Normal, Poisson, Exponential, Gamma, and Uniform distributions. If others are needed, their definitions will given in the quiz. You should also know the definition of moment-generating functions, the characterization of exponential families, and the definition and derivation of the score and Fisher information identities.
- Consider a statistical model given by a family of densities \(f(x|\theta)\) for an outcome \(x\) with parameter \(\theta\) in a open interval in \(\mathbb{R}\). Derive, under standard regularity assumptions, the expected value of the gradient of the log-likelihood evaluated at the ‘true’ value of \(\theta\).
- Consider a statistical model given by a family of densities \(f(x|\theta)\) for an outcome \(x\) with parameter \(\theta\) in a open interval in \(\mathbb{R}\). Derive, under standard regularity assumptions, the identity relating the variance of the gradient of the log-likelihood evaluated at the ‘true’ value of \(\theta\) with the Fisher information.
- Consider the problem of finding an extremum of the real function \(\frac{\log(x)}{1+x} - x^2\). Suppose you are applying the Newton method starting at \(x=1\). Find the next value in the sequence.
- Let the random vector \(x =
\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\) have the following
density for unknown real parameters \(\theta_1, \theta_2 >0\): \[
f(x_1, x_2) = \frac{1}{\theta_1 \theta_2}\exp(-x_1/\theta_1 -
x_2/\theta_2) \quad \textrm{for } x_1, x_2 > 0
\]
- Derive the score function as a function of \(x_1, x_2, \theta_1, \theta_2\).
- Derive the Fisher Information.
- Derive the score function as a function of \(x_1, x_2, \theta_1, \theta_2\).
- Show that the Exponential(\(\lambda\)) distribution for unknown \(\lambda > 0\) belongs to an exponential family but the Uniform(\(0,\theta\)) distribution for unknown \(\theta \in \mathbb{R}\) does not.
Day 17: Monday, November 13
Links:
Sample Quiz Questions for Quiz 4:
Let \(f\) be a twice continuously differentiable real function on \(\mathbb{R}^p\), with a negative definite Hessian. Derive the form of the Newton update to find a maximum of \(f\), given that a Newton update is the maximum of a local quadratic approximation of \(f\).
Consider the problem of finding the root of a real function \(f\) that is twice continuously differentiable on a closed interval \([a,b] \subset \mathbb{R}\) with \(a < b\). Suppose \(f\) has a single root at \(x^* \in (a,b)\). Let \(G(x) = f(x) + x\). Show that \(x^*\) is a fixed point of \(G\) but fixed-point iteration is not a good strategy for finding \(x^*\) if \(|G'(x^*)| > 1\).
Discuss why Fisher Scoring might be sometimes be a better strategy to maximize a log-likelihood than Newton’s method.
Let \(f(x|\theta)\) be a multivariate density for \(x \in \mathbb{R}^p\), with respect to a supporting measure, with the form: \[ f(x|\theta) = \exp \left\{ \sum_{i=1}^p \theta_i x_i - A(\theta) + c(x) \right\} \] with \(\theta \in \Omega \subset \mathbb{R}^p\), where \(\Omega\) is an open set in which the density integrates to 1.
Prove that \(\operatorname{E}(X|\theta) = A'(\theta)\).
Prove that \(\operatorname{Var}(X|\theta) = A''(\theta)\).
Let \(Y\) have a distribution belonging to a scaled exponential family with density of the form: \[ f(y|\theta) = \exp \left\{ \frac{y \theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\} \] with respect to a supporting measure. Let \(\mu(\theta) = E(Y|\theta)\) and let \(\theta\) be modelled as a linear function of two parameters, \(\beta_1\) and \(\beta_2\) with: \[ \theta = z_1 \beta_1 + z_2 \beta_2 \] Show that the gradient of the log-likelihood as a function of \(\beta_1, \beta_2\) has the form \[ \ell'(\beta_1,\beta_2) = \begin{pmatrix} z_1 & z_2 \end{pmatrix} (y - \mu) \]
Day 20: Wednesday, November 22
Links:
Sample Quiz Questions for Monday, November 27
Suppose \((x_i,y_i)\) are pairs of random observations where \(x_i \sim \operatorname{Poisson}(\lambda)\) and \(y_i \sim \operatorname{Poisson}(\alpha\lambda)\), with \(\lambda\) and \(\alpha\) unknown parameters. You have obtained data on 3 pairs, \(i = 1, 2, 3\) but \(y_3\) has gone missing.
Derive the E step and the M step of the EM algorithm to estimate \(\lambda\) and \(\alpha\) with maximum likelihood.
Suppose \((x_i,y_i), i = 1,...,n\) are pairs of random observations where \(x_i \sim \operatorname{Poisson}(\lambda)\) and \(y_i \sim \operatorname{Poisson}(\alpha\lambda)\), with \(\lambda\) and \(\alpha\) unknown parameters. Because of the way the data are obtained, \(y_i\) is missing whenever \(x_i > 3\).
Derive the E step and the M step of the EM algorithm to estimate \(\lambda\) and \(\alpha\) with maximum likelihood.
Do the two problems above when the observations are from the exponential distribution instead of Poisson.
Day 21: Monday, November 27
Day 22: Wednesday, November 29
Day 23: Wednesday, December 4
Sample Questions (largely from the 2nd half of the course):
Consider the following sample test questions in addition to quiz question, midterm test questions and sample questions for the midterm and the quizzes.
Review DAGs and the use of added-variable plots with the Frisch-Waugh-Lovell Theorem to consider the relative values of \(SE(\hat{\beta}_X)\) with different models.
The final exam will cover material from the entire course.
More theoretical:
Suppose you have access to a uniform random number generator. Describe how you could use the Metropolis-Hastings algorithm to generate a Markov chain with a Normal distribution with mean 0 and unit standard deviation as a stationary distribution. Write R code for a function that generates the next observation from this chain.
Let \(h\) on \(\mathbb{R}\) be proportional to a probability density function \(f\), i.e. \(h(x) = c f(x)\) with \(c\) unknown. Let \(g\) be a probability density on \(\mathbb{R}\) that is symmetric around 0, from which you can generate random observations. Describe how you would use random observations with pdf \(g\) to generate a Markov chain with \(f\) as a stationary distribution.
Prove that the process in the previous question has the distribution given by \(f\) as stationary distribution.
Consider estimating the average lifetime of some electronic component, a ‘gizmo’, whose lifetime has an exponential(\(\theta\)) distribution: \[f(x|\theta) = e^{-x/\theta} /\theta \quad x > 0, \theta > 0\] Suppose \(N\) gizmos are observed for a year and \(n\) (\(1 <n < N\)) burn out during the year, with lifetimes recorded as: \(x_1,x_2,...,x_n\) where \(x\) is measured in years (fractional, of course). The remaining \(N-n\) are still functioning at the end of the year.
- Set up the E step of the EM algorithm to estimate \(\theta\). Recall that exponential lifetimes are ‘memoryless’, i.e. the distribution of the additional life of a gizmo that survives one year is exponential(\(\theta\)).
- Derive the M step, i.e. how do you obtain \(\theta_{t+1}\) as a function of \(\theta_t\)?
Prove that each step in the EM algorithm cannot reduce the log-likelihood. Assume that the densities in the model, \(f(x|\theta)\) are positive with the same support.
Assume known as a general fact that the Kullback-Leibler divergence from the distribution \(f(\cdot|\psi_1)\) to \(f(\cdot|\psi_2)\) is non-negative, i.e. \[\operatorname{E}\{\log F(x|\psi1_1) - \log F(x|\psi_2) | \psi_1\} \ge 0\]Describe the Newton-Raphson method for finding the MLE of a statistical model. Give a sketch of a proof, using Taylor series expansions, that it converges to the maximum if the starting value is sufficiently close to the maximum.
Describe the Newton-Raphson and the Fisher Scoring method for finding MLEs and compare some advantages and disadvantages of each method.
Describe the Nelder-Mead algorithm in detail. (You’ll forget the details after the exam but it’s worth having to learn it at some point in your career.) Be ready to apply it for one or two steps in a simple example.
Applied:
The following 4 questions are based on the ‘Vocab’ data set used in
Assignment 3.
Recall that it records a vocabulary score for over 30,000 subjects
tested over the years between 1974 and 2016. The questions below refer
to the following output in R. Note that the variable ‘Year’ is measured
in years since 2000.
library(car)
summary(Vocab) year sex education vocabulary
Min. :1974 Female:17148 Min. : 0.00 Min. : 0.000
1st Qu.:1987 Male :13203 1st Qu.:12.00 1st Qu.: 5.000
Median :1994 Median :12.00 Median : 6.000
Mean :1995 Mean :13.03 Mean : 6.004
3rd Qu.:2006 3rd Qu.:15.00 3rd Qu.: 7.000
Max. :2016 Max. :20.00 Max. :10.000data <- within(Vocab,
{
Year <- year - 2000
Year2 <- Year^2
Year3 <- Year^3
Yearf <- factor(Year)
}
)
fit1 <- lm(vocabulary ~ factor(Year) * sex, data = data)
fit2 <- lm(vocabulary ~ (Year + Year2 + Year3) * sex, data = data)
fit3 <- lm(vocabulary ~ Year * sex, data = data)
summary(fit2)
Call:
lm(formula = vocabulary ~ (Year + Year2 + Year3) * sex, data = data)
Residuals:
Min 1Q Median 3Q Max
-6.1730 -1.0826 0.0162 1.0994 4.1273
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.141e+00 2.790e-02 220.123 < 2e-16 ***
Year 1.215e-02 2.793e-03 4.351 1.36e-05 ***
Year2 -8.358e-04 1.880e-04 -4.445 8.83e-06 ***
Year3 -4.831e-05 1.017e-05 -4.749 2.06e-06 ***
sexMale -1.203e-01 4.226e-02 -2.847 0.00441 **
Year:sexMale -4.446e-03 4.230e-03 -1.051 0.29331
Year2:sexMale 3.865e-04 2.842e-04 1.360 0.17376
Year3:sexMale 2.488e-05 1.533e-05 1.623 0.10461
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.114 on 30343 degrees of freedom
Multiple R-squared: 0.001568, Adjusted R-squared: 0.001337
F-statistic: 6.806 on 7 and 30343 DF, p-value: 4.252e-08summary(fit3)
Call:
lm(formula = vocabulary ~ Year * sex, data = data)
Residuals:
Min 1Q Median 3Q Max
-6.0848 -1.0590 -0.0145 1.0799 4.1114
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.0504270 0.0172503 350.742 < 2e-16 ***
Year 0.0021493 0.0012801 1.679 0.09315 .
sexMale -0.0799758 0.0260553 -3.069 0.00215 **
Year:sexMale 0.0009994 0.0019247 0.519 0.60359
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.115 on 30347 degrees of freedom
Multiple R-squared: 0.0006394, Adjusted R-squared: 0.0005406
F-statistic: 6.472 on 3 and 30347 DF, p-value: 0.000225anova(fit3, fit2, fit1) Analysis of Variance Table
Model 1: vocabulary ~ Year * sex
Model 2: vocabulary ~ (Year + Year2 + Year3) * sex
Model 3: vocabulary ~ factor(Year) * sex
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30347 135692
2 30343 135566 4 126.03 7.0624 1.116e-05 ***
3 30307 135208 36 358.24 2.2305 3.234e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pred <- within(
with(data, expand.grid(Year = unique(Year), sex = unique(sex))),
{
Year2 <- Year ^ 2
Year3 <- Year ^ 3
}
)
pred$fit1 <- predict(fit1, newdata = pred)
pred$fit2 <- predict(fit2, newdata = pred)
pred$fit3 <- predict(fit3, newdata = pred)
library(latticeExtra) Loading required package: latticexyplot(fit1+fit2+fit3 ~ Year, pred, groups = sex, type = 'b')
Is the test whose p-value appears in the third line of the output for the
anovafunction, a legitimate hypothesis test. Why or why not? Describe in words what the p-value is indicating, if anything.Using the model
fit2, estimate the rate at which the gender gap is estimated to be changing in the year 2010. How would you go about obtaining a confidence interval for this estimate?Answer the previous question using the model
fit3. Which model do you think provides a better answer to this question? Why?Describe the interpretation of the estimated coefficients for
sexMaleand forYear:sexMalein the summary output forfit3.Describe how you would simulate an observation from a distribution with density \(f(x) = exp\{-x/2\}/2, x >0\) using the inverse CDF. Write the code for a function that would generate \(n\) observations from this distribution using
runifas a uniform random number generator.
References
R is to S as Linux is to Unix as GNU C is to C as GNU C++ is to C …. S, Unix, C and C++ were created at Bell Labs in the late 60s and the 70s. R, Linux, GNU C and GNU C++ are public license re-engineerings of the proprietary S, Unix, C and C++ respectively.↩︎