Day 1: Monday, January 6

• Lecture slides and Quiz 1 questions

• Course description

• This evening, I will resend invitations to join Piazza to all students currently registered in the course.

  • I will use your York e-mail address. If you don’t read email at your York email address, please make sure that it’s forwarded to an email address that you do read regularly.
  • Please do not change your e-mail address on Piazza because your York email address is used to identify you so you get credited for your contributions. You make my life much easier if I don’t have to manually cut and paste your grades.

• Topic 1: AMSTAT News Survey of Master’s Graduates (2023 US Statistics Graduates surveyed in mid-2024)

  • 

• Topic 2: The meaning of p-values Consider the simplest of statistical tests. Suppose you plan to sample 100 observations from a population that is known to have a normal distribution with variance equal to 1 but unknown non-negative mean $\mu$.

  You plan to test $H_0: \mu = 0$ versus the one-sided alternative $H_A: \mu > 0$.

  Consider what happens if you use a 5% test, i.e. reject the null hypothesis if $p < 0.05$.

  Suppose you reject $H_0$. What does this tell you about $H_0$?

  What if the $p$-value is 0.049? Or what if the $p$-value is 0.0003?

Day 2: Wednesday, January 8

• Topic 2: How harmful is X?
  • The Trapezoid of Means: Simpson’s Paradox, the Paik-Agresti diagram and the Liu-Meng diagram:
    • Is using Sun Screen Lotion bad for you? R script
      • Try the exercise at the end of this script.
  • (On Day 3) How harmful is smoking cigarettes?: Here’s a link to the R script

Day 3: Friday, January 10

Announcements:

• Tutorial today at 4:30 on Zoom. See link under Quick Links

• Tutorials on Tuesday mornings at 9am on Zoom. See link under Quick Links

• See modified course description

• Topic 1: The Trapezoid of Means: Simpson’s Paradox with Paik-Agresti and Liu-Meng Diagrams / R script

• Topic 2: How to estimate everything

• Topic 3: Regression Review: Regression in R

Day 4: Monday, January 13

Quiz on Wednesday:

Sample question:

Given a table of conditional mean survival rates in percents(Y) by Gender (G: F or M) for two treatments (X: A or B) for a disease D:

The following are the numbers of individuals in each group:

1. Draw the trapezoid of mean survival rates.
2. Find the marginal effect of X (comparing treatment B to treatment A)
3. Find the conditional (specific) effects of X
4. Find the marginal effect of G.
5. Draw the lines corresponding to the Paik-Agresti diagram and the Liu-Meng diagram.
6. Would the conditional effects or the marginal effect be more appropriate as a measure of the effectiveness of this treatment assuming the Ms and Fs were assigned randomly (albeit not with equal probability) to the treatments.
7. Suppose variable G stands for ‘Gastric Complications’ that can be caused by the treatment (F is few and M is many). Again, subjects (a random mix of men and women) were randomly assigned to the the two treatments. Would the conditional effects or the marginal effect be more appropriate as a measure of the effectiveness of this treatment.
8. Discuss your answers to (f) and (g) briefly.

Day 5: Wednesday, January 15

Announcement: Tutorial on Zoom (click here), every Tuesday from 1:30 to 2:30 pm.

Quiz Today:

• Solution

• Topic 1 (continuing): Regression Review: Regression in R / R script

• Topic 2: Some real data on Smoking

• Topic 3: Simple Regression to Multiple Regression

Day 6: Friday, January 17

Concepts and Theory

• Topic 1: (continuing) Simple Regression to Multiple Regression
• Topic 2: Multiple Regression in $\beta$-space
  
• Topic 3: Three Basic Theorems
• Topic 4: Causal Questions
• Topic 5: The Causal Zoo
• Topic 6: 23 statements about Statistics (one is correct)

Using R:

• Regression Review: Regression in R / R script
• Working with Data in R / R script

Learning R:

Why R? What about SAS, SPSS, Python, 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, CRAN had 22,243 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, ² 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 1/2 million 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.

Copy the following scripts as files in RStudio:

• Elementary examples in R: CAR_1.R, and
• More advanced techniques: Working with Data.R.

Play with them line by line.

Post questions arising from these scripts to the ‘question’ folder on Piazza. We will take up some questions in class and others in tutorials scheduled to deal with questions on R.

Day 7: Monday, January 20

Day 8: Wednesday, January 22

• No class this day because of an unfortunate conjunction of unlikely events.

Day 9: Friday, January 24

Day 10: Monday, January 27

Sample question for Quiz 3 on Wednesday, January 29:

• Sample questions
• See also the question below.

Topic 1:

• Recap so far
• What plots do the work of simple scatterplots for multiple regression? Let your imagination run wild depending on the nature of the process generating the data, the shape of the data, and the purpose of the analysis. Here are two basic plots:

  1. Added-Variable-Plot (AVP):
     Plot $Y - \hat{Y}|_{other X's}$ against $X - \hat{X}|_{other X's}$

     • Shows the first-order leverage and influence of each observation on $\hat{\beta}_i$

     • Note that the definition of the AVP yields the simple scatterplot centered at the origin if the only ‘other variable’ is the intercept.

     • In R use:

       library(car)
fit <- lm(income ~ education + prestige, data = Prestige)
avPlots(fit, ellipse = TRUE)
?avPlots # for more options

     When you look at this plot, think of how you would interpret it if it were a simple scatterplot.
     You can actually construct a 95% confidence interval for the multiple regression coefficient, $\hat{\beta}_i$ for the effect of $X_i$ keeping other $X$’s constant, using this plot the same way you used a simple scatterplot. Of course, the 95% confidence interval also gives you the result of a 5% test of the hypothesis that $\beta_i = 0$ in the multiple regression.

     Note: According to Wikipedia, the person who first discovered this fact is a statistician, Udny Yule, who published it in 1911. Nevertheless the theorem is known as the Frisch-Waugh-Lovell Theorem. The eponymous economist, Ragnar Frisch, won the first Nobel Prize in Economics in 1969. Udny Yule is also considered the ‘real’ discoverer of Simpson’s Paradox, many years before Simpson published his paper on it. Thank goodness theorems and paradoxes are rarely named after their original discoverers. E.g. we might have no idea what someone is talking about when they refer to ‘Yule’s Theorem’. I hope he considers that a consolation.

  2. Residual versus Leverage plots:
     Plot the standardized residual $\frac{e_i}{\hat{\sigma} \sqrt{1 - h_{ii}}}$ against $h_{ii}$.
     Shows the first-order leverage and influence of each observation on $\hat{Y}_i$.
     In R use:

     plot(rstudent(fit) ~ hatvalues(fit))

     or look at the fourth plot when you just use

     plot(fit)

     This plot will also show you selected isoquants of Cook’s Distance, one of perhaps over 50 diagnostics generated by two groups in the late 70s: David Cook and Sanford Weisberg in one group, and David Belsley, Edwin Kuh and Roy Welsch, in the other. Some analysts look at tables of these diagnostics and apply conventional cutoffs to flag problematic points. Many diagnostics are functions of each other and I think that we get much better insights by looking at them graphically, allowing us to detect isolated outliers or clusters of outliers seeing patterns involving more than one numeric diagnostic simultaneously.

     Note that $$h_{ii} = \frac{1}{n}\left( 1 + \Delta_i^2 \right)$$ where $\Delta_i$ is the Mahalanobis distance of the i-th point from the predictor values in the sample. In simple regression, this is just the ‘z-score’ for $X$ for the i-th observation.
     $h_{ii}$ obeys the following inequalities if the model has an intercept: $$\frac{1}{n} \le h_{ii} \le 1$$ and $h_{ii} = 1$ iff all the points omitting the i-th point are perfectly collinear. What inequalities does this imply for $\Delta_i$? What is the implication if $h_{ii} = 1/n$?
     Quiz question: Where would you expect to find outliers of ‘Type 1’, ‘Type 2’, and ‘Type 3’ in this plot.

Topic 2 (continuing): Regression Review

• Simple to Multiple Regression
  • Pay equity in a small law firm: Simpson’s Paradox, Paradox lost
    
• Multiple Regression in $\beta$-space
  
• Three Basic Theorems

Day 11: Wednesday, January 29

Day 12: Friday, January 31

Day 13: Monday, February 3

Day 14: Wednesday, February 5

The class today was on the blackboard. See the recording for its content.

Day 15: Friday, February 7

• Sample questions for the midterm test

• Quiz 3 solutions

• Summary on causality using linear causal DAGS

  • Backdoor criterion to estimate the causal effect of X on Y: (See also the last page of the link above)
    If your linear causal graph and model satisfy both of the following conditions, then you have satisfied sufficient conditions so that your estimate of the effect of X is an estimate of a causal effect of X on Y:
    1. Trace out all backdoor paths between X and Y and ensure that all backdoor paths are blocked. A backdoor path is blocked if either/or:
       1. There is a collider on the path that is NOT controlled. ‘Controlled’ can mean that the variable is included in a regression model, or that the analysis is stratified on the variable, or – more subtly and challenging to detect – that the data was selected on the basis of the variable and only one or a few levels are observed, resulting in selection bias.
       2. One or more variables that are not colliders on the path are controlled. If a path has been blocked by controlling for a collider variable, it can be unblocked by controlling for other non-colliders on the path.
    2. Do not control any descendants of X. That includes, of course, descendants of Y.

  So, step-by-step:

  1. Trace out the backdoor paths.
  2. Make sure all the backdoor paths are blocked using the criteria above: either a collider that is not controlled or one or more non-colliders that are controlled. If a collider is controlled (through the analysis by inclusion in the model, or through selection) then the path, thus unblocked, can be re-blocked by controlling one or more non-colliders.
  3. Make sure the model does not include descendants of X.

  Note that average causal effects can only be generalized to members of the population that were included in the observed sample. For example, if the sample was selected through a collider variable and the resulting backdoor path re-blocked by the inclusion of another variable along the path, then the observed average causal effect can only be generalized to members of the population that were selected through the collider.

• Recommended Piazza weekly feedbacks:

  • An extended discussion of colliders
  • What if we can’t come up with a reasonable DAG?
  • A DAG in Wednesday’s lecture on the blackboard
  • What’s accomplished by a dummy variable for an outlier?

• To further explore DAGs and related concepts:

  • Download the dagitty package in R and play with it. Visit the home page for the project. A DAG implies certain relationships among its variables. Some of these relationships are purely causal and can’t be checked with a statistical test on the data. For example, with the simple DAG we used for Coffee, Stress and Heart Damage, the correct causal analysis depended on the direction of the arrow between Coffee and Stress. If the arrow points from Coffee to Stress, then Stress is a mediating variable and should be excluded from the regression. If the arrow points from Stress to Coffee then Stress is a confounding factor and should be included in the regression.
    There is no statistical test to check which way the arrow should go because either direction gives rise to the same means, variances and covariances among the three variables, and that’s all a statistical test would see.
    However, with more complex DAGs, the DAG may have both causal and correlational implications. You can use the dagitty package to test whether the correlational assumptions are satisfied and, if they’re not, you can explore further to refine your DAG.
    Just playing with examples from the dagitty package should be both fun and very instructive in deepening your understanding of causality.

• Recommended books to go deeper into causal models:

  • Hernan and Robins (2020) Causal Inference: What If? (hernan2020?)
  • Pearl and Mackenzie (2018) The Book of Why (Pearl and Mackenzie 2018)

Day 16: Monday, February 10

Review

Day 17: Wednesday, February 12

Midterm Test

Day 18: Friday, February 14

Mixed Models:

• Hierarchical Models
• R script for hierarchical models: Lab 1.R
• Longitudinal Models
• R script for longitudinal models: Lab 2.R

Day 19: Monday, February 24

Class links:

• Small hierarchical example comparing pooled, between, within and mixed models, with and without a contextual variable
  • R script for the above
  • Means only model: Mixed Models or Pooled Analysis (rough draft)
  • When should you not use a contextual variable for a lower level predictor:
    • When the predictor is ‘balanced’, e.g. has the same values, or the same mean in each cluster: In that case the contextual variable may be almost constant and collinear with the intercept.
    • When it’s reasonable to believe that the between-effect is equal to the within-effect – and there is no evidence to the contrary – so it’s reasonable to ‘borrow strength’ from the between effect in the estimation of the within-effect.
• Multilevel Mixed Models
  • R script to play with Multilevel Models: Lab 1.R
• Non-Linear Mixed Models
• Longitudinal Linear Models
  • R script to play with Longitudinal Models: Lab 2.R
    • This is a good examples of why models matter and how the estimated comparisons between drugs depend on which potential confounding factors are included in the model. This includes longitudinal components such as trends and possible autocorrelation between measurements that are close to each other in time. As you had more components the story keeps changing.
      
• Added-Variable Plot (Frisch-Waugh-Lovell Theorem) and the Linear Propensity Score Theorem
  • Provides important insights into equivalent models and nearly equivalent models
• Parametric splines
  • Models that are non-linear in X but linear in parameters.
• Dealing with Heteroskedasticity

Day 20: Wednesday, February 26

Day 24: Friday, March 7

• Guide to Parametric Splines
• Longitudinal Parametric Splines: Sleep at Birth sample script

…

Day 27: Friday, March 14

• Non-Linear Hierarchical Models
• Lab for non-linear mixed models

…

Day 30: Friday, March 21

• Revisiting the G matrix and the Role Centering Predictors within Groups (CWG)
  • Longitudinal Models
  • Contextual and Compositional Effects

Day 33: Friday, March 28

• Longitudinal Models
• Principle of Marginality
• Sample final exam questions
  • This sample final has many more questions than the actual final. The final exam, lasting 2 hours, will have questions similar to these but whose marks sum to 100.

Day 34: Monday, March 31

• 23 Statistical Statements
  • Which ONE is TRUE?
• When, why and how to use REML or ML when fitting mixed model?
  • REML vs ML
  • REML vs ML R script