Disclaimer

• I am not a statistician and do not pretend to fully understand all the implications of what will be presented.
• Some sloppy mathematical notation is used. If you know how to improve it, please suggest some modifications.
• Not all the code is shown in the presentation to reproduce the examples.
  • You have to go in the
  • presentation RMarkdown files or
  • the script file Workshop R script, which requires the 0.0_initialize.R and other scripts specified within.
  • Don't forget to change the paths associated with some files to load everything.
• In addition, this is a work in progress. Any input you have is more than welcome.

Outline

Theory:

1. Why are simulations relevant?
2. When simulations are useful?
3. What you need to know before doing simulations?
4. Let's explore what underlie some random processes
5. Statistical distributions in more details
6. Simulating data for models (lm, glm, lmm, in development)
7. Simulating data for power analysis (in development)
8. Advanced simulations
9. References
10. Acknowledgements

The exercices for this workshop: Challenge wallpaper, Challenge 1, challenge 2, challenge 3, challenge 4, challenge 5, challenge 6, challenge 7.

Supporting material:

Please check out the Simulation in R Cheat Sheet

Prerequisites

For this workshop it is useful to have a solid understanding of

1. Linear models and generalized linear models (GLMs)
2. Programming in R

See the QCBS R workshops if you want to revise these topics. Select the R workshops needed.

• We are going to review some aspects of these workshops, but getting some experience from those subjects can clearly help you for this workshop.
• We are going to be more explicitly show assumptions of models and use for loops in this workshop
• You should also have installed R on CRAN here and you could get RStudio as well.

Learning objectives

1. Develop intuition and skills to perform simulations
2. Use specific and relevant functions in R when designing simulations
3. Articulate what are the assumptions underlying a simulation (for the models tested or the scope of the simulation)
4. Know how and when to use important statistical distributions
5. Simulate different data structure to perform LMs, GLMs, and (G)LMMs.

Expectation from the course

• Please, I invite you to participate as much as you can.
• I certainly have a lot to learn and I wish that our learning environment is welcoming to committing errors and learning from them

Background and questions

• Have you ever performed a simulation in R?
• What are simulations for? Would simulations be useful for your current or future work?
• Why do you think they are important?
• When do you think they are helpful?
• What is a model? Can you give examples of models that are not statistical models?
  

Background and questions

• Biological models

  • Exponential growth: $n(t+1)=Rn(t)$ (discrete) or $\frac{dn}{dt}=rn(t)$ (continuous)
  • Logistic growth: $n(t+1)=n(t)+rn(t) \Big(1-\frac{n(t)}{k} \Big)$ (discrete) or $\frac{dn}{dt}=rn(t) \Big(1-\frac{n(t)}{k} \Big)$ (continuous)
  • Haploid and diploid selection
  • Competition equations
  • Consumer-resource equations
  • SIR equations (Susceptible, Infected, Recovered)

• Statistical models

  • LM, GLM, GLMMs, GAM, Ordinations, etc.
  • Machine learning
  • Dynamical models

Survey: defining models

Usefulness of models

"Models can help to guide a scientist's intuition about

• how various processes interact;
• they can point out logical flaws in an argument;
• they can identify testable hypotheses, generate key predictions, and suggest appropriate experiments; and
• they can reshape fields by providing new ways of thinking about a problem.

But mathematical models also have their limitations. The results of an analysis are only as interesting as the biological questions motivating a model. [...]"

— Sarah. P. Otto, and Troy Day 2007 (modified)

Survey: defining simulations

Defining simulations

— Byron J. T. Morgan 1984

Survey: would simulations be useful for you?

Simulations are games: we make up the rules

• A challenges is to figure out what parameters (or rules) are governing a natural process.
• Simulations are powerful tool that can generate 'alternative worlds' or models of the world.
• We are responsible and need to be conscious about the choices we make in order to get an answer.
• In summary:
  • We try to understand natural processes, which are often manifestations of random variables or stochastic process
• Being able to do simulations can be a nice addition to your research toolkit.

Are these two distribution showing the same stochastic process?

Why simulations? Find patterns from stochastic processes

• Simulations are useful to test the properties of randomly generated data
• Since we design the simulation, we know the parameters of the processes.
• It is then possible to test various methods or statistical analysis to
  1. learn how data are generated,
  2. Understand and see a model works and verify its assumptions or predictions,
  3. explore old and new ideas, theories, hypothesis, the existence of a phenomenon or a slow process,
  4. prepare your analysis script before data collection (or preregistration),
  5. do power analysis and dry run an experimental design,
  6. simulate tragic events to see their impact (e.g., nuclear explosion, flood, emergence and trajectories of hurricanes, etc.)
  7. give you teaching material,
  8. etc.

Pros and cons of simulations

The list below was taken from Banks (1998 Handbook of simulation).

ALL THE TIME!!!!

• Simulations can help you at all stages of your research from planning to publishing and beyond!
• Whether you want to:

1. Learn about science in general

   • learn how a mechanism work (i.e., genetic drift, natural selection),
   • develop a deeper understanding of a concept.

2. Design a particular study

   • Specify and refine a research objective, question, hypotheses to test and their predictions,
   • a research plan, an experimental and sampling design or test a statistical method,
   • play with data before your field work,
   • simulate random point on a map to sample it.

3. Publish an article or criticize

   • respond to reviewers for a paper you want to publish,
   • support an argument or a point you want to make,
   • verify the claims made by other researchers.

4. etc.

Quick example: Genetic drift without mutation

• Here we have 4 plots showing genetic drift with the same simulation.
• The only difference is the number of individuals in each population
• Your can deduce many things from the allele frequency change.

Quick (?) example: Natural selection

• How would you model natural selection?
• What is the type of selection we want to simulate?
• What is being selected? What is that type of data (continuous [range?], discrete [counts?])?

• But what all of these mean?
• What do they assume?
• What is the data generated to make these lines and curves?

What do you see here?

• Can you give as much information as you can within and between each graphs.

Find correlation and covariance

• The correlation of the first graph is 0.92 and the covariance is 45.74.
• The correlation of the second graph is 0.97 and the covariance is 2.04.

What do you see here?

• Can you give as much information as you can within and between each graphs.

• Here the mean of X and Y is shown

What do you see here?

• Can you give as much information as you can within and between each graphs.

• Here the mean and the standard deviation of X and Y is shown

What do you see here?

• Can you give as much information as you can within and between each graphs.

• Here the mean and the standard deviation of X and Y is shown with the regression lines

What do you see here?

• Can you give as much information as you can within and between each graphs.

• Here the mean and the standard deviation of X and Y is shown with the regression lines and the distributions of X and Y.

What do you see here? (2)

• Can you give as much information as you can within and between each graphs.

What do you see here? (2)

• Can you give as much information as you can within and between each graphs.

What do you see here? (2)

• Can you give as much information as you can within and between each graphs.

Tips when performing simulation

Description section. Add a 'Description' section to all your simulation scripts to introduce what the script is about.

• You can add information how to use important arguments or a step-by-step description on how to use the script.
• You can add some references that you used to build your scrip.

# Description  ------------------------------------------------------------
#### ### ### ## #### ### ### ## #### ### ### ## 
# TITLE OF THE SCRIPT
# Created by YOUR_NAME
# 
# Why: 
# Requires:
# NOTES: 
# Reference : 
#### ### ### ## #### ### ### ## #### ### ### ## 
# Code here ...

Tips when performing simulation

Comments. Be extra generous when commenting your code to describe as precisely as possible for the information that is not variable.

# Description  ------------------------------------------------------------
# This is the description section as previously presented 
# Libraries ---------------------------------------------------------------
# Here load the libraries used in the script 
library(ggplot2)
# Functions ---------------------------------------------------------------
# Add and describe the functions used in the script 
## Add more information
function.name = function(var){tmp=2+var; return(tmp)}
# Plots -------------------------------------------------------------------
# Plotting the data simulated 
## Add more information
plot(function.name(1:10)+rnorm(10))

(We used fewer comments to lighten the presentation, but please add comments.)

Tips when performing simulation

• Structure your work space in a way that is easy to track:
  1. README.txt # explain the project and minimal details on how to navigate the project, and metadata
  2. data # Place all the data in this folder
  3. raw # it is suggested to 'lock' the raw data, so you don't accidentaly modify the data
  4. processed # This would be the 'clean' data
  5. output # This is what you generate in R
  6. plots # Export the plots with pdf, png, svg, or anyother function to save your plots
  7. data # This would be intermediate data, tables, and other intermediate manipulation of data
  8. scripts # In this folder, place the scripts (You can have a 'initialize.R' script and 'runall.R' script to be able to reproduce the analysis. Check out the sessionInfo() function)
  9. functions # You might want to generate your own custom functions. They can be placed here or in the 'initialize.R' script.
  10. temp # If needed, place temporary files and want to delete them after

Building simulations is an iterative process

ATTENTION. Usually, coding simulation is written in small increments and then put together in a cohesive function that will make what we are looking for.

• It takes practice to code simulations: the more you code, the more you develop tricks.
• Always try to relate your simulation to a certain question of interest. This will bring more meaning to the insights that you are trying to make.

In this section, we will discover functions that might help when coding simulations in R.

Functions useful in simulations (RNG)

• When performing simulations, you will have to play with (pseudo-)randomly generated numbers from a random number generator (RNG).
• This is a challenge if we want to replicate the analysis we are performing.
• R has the function set.seed() that help us to play with the RNG.

The example below uses a RNG to extract numerical value between 1 and 10

runif(n = 1, min = 1, max = 10) # Gives a random number between 1 and 10
# [1] 3.588198
runif(n = 1, min = 1, max = 10) # RNG wasn't reset, different answer (see above)
# [1] 8.094746
runif(n = 1, min = 1, max = 10) # Different again... 
# [1] 4.680792
set.seed(42); runif(n = 1, min = 1, max = 10) # This sets the RNG 
# [1] 9.233254
set.seed(42); runif(n = 1, min = 1, max = 10) # The exact same number 
# [1] 9.233254

Functions useful in simulations (sample)

• The function sample() randomly picks a value from a vector (i.e., random sampling).

The example below uses a RNG to extract numerical value between 1 and 10

set.seed(12) # Set the RNG 
v.1.10 = 1:10 # Make a vector from 1 to 10 
# Randomly pick 1 (size) value from the vector (x), without replacement 
sample(x = v.1.10, size = 1, replace = FALSE) 
# [1] 2

• The values don't have to be numerical: they could be characters or factors

set.seed(3) # Set the RNG 
# Randomly pick 5 (size) letters from the vector (x), without replacement 
sample(x = LETTERS, size = 5, replace = FALSE) 
# [1] "E" "L" "G" "D" "H"
sample(x = as.factor(month.abb), size = 5, replace = FALSE) 
# [1] Nov Aug Apr Jul Dec
# Levels: Apr Aug Dec Feb Jan Jul Jun Mar May Nov Oct Sep

Functions useful in simulations (sample)

• The function sample() can actually be used in order to do permutations
• Let's say we have a data frame with 2 columns

set.seed(123)
n = 40; col = viridis::viridis(n = n)
x = 1:n ; y = 2+.5*x + rnorm(n, sd=7)
df.xy = data.frame(x,y, col )

set.seed(321)
df.xy$x.s = sample(df.xy$x)
df.xy$y.s = sample(df.xy$y)
# We break up the link of X and Y

Functions useful in simulations (sample)

• Below is an example using sample to generate
  • permutation of the data (orange lines) and
  • Bootstrap, leave-'one'-out data points (here leaving 20 out), recalculating the lm (red lines)

permutate.df = replicate(n = 200, # nperm 
          expr = data.frame(y = sample(df.lm$y,size = nrow(df.lm), replace = FALSE),
                            x = df.lm$x), simplify = FALSE)
bootstrap.fun <- function(data) {
  tmp.sample = sample(1:nrow(data),size = c(nrow(data)-20), replace = FALSE)
  data.frame(y = data[tmp.sample,"y"], x = data[tmp.sample,"x"]) } # end boot fun
bootstrap.df = replicate(n = 200, expr = bootstrap.fun(df.lm), simplify = FALSE)

Functions useful in simulations (sample)

• You can also simulate random dates in a time interval.

set.seed(42)
# Get a sequence of dates
datae.seq = seq(from = as.Date('2010/01/01'), 
                to = as.Date('2022/01/01'), 
                by = "day")
# Look at beginning and end of sequence 
head(datae.seq, 4); tail(datae.seq, 4)
# [1] "2010-01-01" "2010-01-02" "2010-01-03" "2010-01-04"
# [1] "2021-12-29" "2021-12-30" "2021-12-31" "2022-01-01"
# Get only 5 elements of the generated sequence 
sample(datae.seq, 5)
# [1] "2017-02-21" "2021-02-20" "2016-06-26" "2013-01-02" "2013-06-05"

Functions useful in simulations (sample)

• You can also generate a die and randomly toss it.
• You can change the probability of landing on one face more than the other.

set.seed(50)
p_dice = c(1,1,1,1,1,5) # Here we have 5 times the change of landing on 6
                        # Same as writing p_dice/sum(p_dice) or the prob.
nb.tosses = 100
die_results <- sample(x = 1:6, # or seq(from = 1, to=6, by=1)
                      size = nb.tosses,
                      replace = T, prob = p_dice) 
barplot(table(die_results), ylab = "Fq", xlab ="Face", main = "Loaded dice") # table(die_results)/nb.tosses

Functions useful in simulations (rep)

• Numerical values are not the only type of class you'll see in simulations. Characters (or factors) can be generated in R easily.
• The rep() function can help you with this

(let4first = LETTERS[1:4])
# [1] "A" "B" "C" "D"
rep(let4first, times = 2) # Repeat the WHOLE sequence twice 
# [1] "A" "B" "C" "D" "A" "B" "C" "D"
rep(let4first, each = 2) # Repeat each element twice 
# [1] "A" "A" "B" "B" "C" "C" "D" "D"
# Set the number of repeat for each element in the vector 
rep(let4first, times = c(1,2,3,6))
#  [1] "A" "B" "B" "C" "C" "C" "D" "D" "D" "D" "D" "D"
# Complete replication: replicate each element twice and do this three times 
rep(let4first, each = 2, times = 3)
#  [1] "A" "A" "B" "B" "C" "C" "D" "D" "A" "A" "B" "B" "C" "C" "D" "D" "A" "A" "B"
# [20] "B" "C" "C" "D" "D"
rep(let4first, length.out = 6) # Repeat the vector until you hit the length.out
# [1] "A" "B" "C" "D" "A" "B"

Functions useful in simulations (gl)

• There are functions in R that can help you generate factor levels (with the gl function).

nb.of.levels = 2
nb.of.replicates = 8
labels.for.the.factors = c("Control", "Treat")
## First control, then treatment:
gl(n = nb.of.levels, k = nb.of.replicates, labels = labels.for.the.factors)
#  [1] Control Control Control Control Control Control Control Control Treat  
# [10] Treat   Treat   Treat   Treat   Treat   Treat   Treat  
# Levels: Control Treat
## 20 alternating As and Bs
gl(n = 2, k = 1, length = 20, labels = LETTERS[1:2])
#  [1] A B A B A B A B A B A B A B A B A B A B
# Levels: A B
## alternating pairs of 1s and 2s
gl(n = 2, k = 2, length = 19) # see last element missing
#  [1] 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2
# Levels: 1 2

Functions useful in simulations (replicate)

• In some cases, you want to replicate data, for example when you want to generate multiple populations at once with the same parameters
• replicate() can be used instead of a for loop and make simulations faster.

set.seed(1234)
data.replicated = replicate(n = 2,
                            expr = data.frame(gr = rep(LETTERS[1:3], each = 2),
                                              y = rnorm(6)), 
                            simplify = FALSE)

data.replicated[[1]]
#   gr          y
# 1  A -1.2070657
# 2  A  0.2774292
# 3  B  1.0844412
# 4  B -2.3456977
# 5  C  0.4291247
# 6  C  0.5060559

data.replicated[[2]]
#   gr          y
# 1  A -0.5747400
# 2  A -0.5466319
# 3  B -0.5644520
# 4  B -0.8900378
# 5  C -0.4771927
# 6  C -0.9983864

Functions useful in simulations (expand.grid)

• The function expand.grid() will make a data frame from all combinations of factor variables
• It can be useful when making hypothesis testing or want to explore an experimental design

set.seed(1234)
exp.design = expand.grid(height = seq(60, 70, 10), 
                         weight = seq(100, 200, 100),
                         sex = c("Male","Female"), 
                         stringsAsFactors = TRUE) # characters will be factors
exp.design
#   height weight    sex
# 1     60    100   Male
# 2     70    100   Male
# 3     60    200   Male
# 4     70    200   Male
# 5     60    100 Female
# 6     70    100 Female
# 7     60    200 Female
# 8     70    200 Female

Functions useful in simulations (expand.grid)

• For example, you could expand.grid() to make a deck of cards

cards = c( "A",2:10, "J", "Q", "K"); length(cards) # Get cards type
# [1] 13
suits = c("♦", "♣", "♥", "♠" ) # Get the suits c("Diamonds", "Clubs", "Hearts", "Spades")
cute.deck <- expand.grid(cards = cards, suits = suits)
cute.deck$nice.cards = apply(cute.deck, 1, function(x) paste0(x, collapse = "")) # Combine cards and suits
cute.deck$col = ifelse(cute.deck$suits %in% c("♦","♥"),"red","black") # add colour 
# Select cards at random 
set.seed(1234)
n.c = 5
row.sel = sample(1:nrow(cute.deck),size = n.c, replace = FALSE)
cute.deck[row.sel,"nice.cards"]
# [1] "2♥" "3♣" "9♣" "J♥" "5♠"

Functions useful in simulations (outer)

• The function outer() can apply a function to the input arrays.

x <- 1:10; names(x) <- x
# Multiplication & Power Tables
x %o% x # same as outer(x, x, "*")
x.mul = outer(x, x, "*")
x.sub = outer(x, x, "-")
x.add = outer(x, x, "+")
x.div = outer(x, x, "/")

Challenge 0: wallpapeR

• Provide 2 functions (see the #<<) with their argument(s) to

1. store a random number between 1 and 115 in rdm.nb and
2. one that will give the same number from the pseudo-random generator each time the script is ran.

#<<
par(bg = NA, mar = c(0,0,0,0)) # will take the whole screen
x = 1:70; y = 1:70 # make variables 
x.val = (x^2);   y.val = y
xy.out = outer(x.val, y.val, "-") # manipulate the variables
nb.col.pal = length(grDevices::hcl.pals()) # find all palettes in hcl.pals
#<<
xy.out = xy.out + rnorm(length(xy.out), mean = 0, sd = 200) # Prettify with RDM
image(t(xy.out),xaxt = "n", yaxt = "n", bty = "n",
      col = hcl.colors(1000, palette = hcl.pals()[rdm.nb])) # 61,93 looks good

Challenge 0: wallpapeR - Solution

• Provide 2 functions (see the #<<) with their argument(s) to

1. store a random number between 1 and 115 in rdm.nb and
2. one that will give the same number from the pseudo-random generator each time the script is ran.

set.seed(9)
par(bg = NA, mar = c(0,0,0,0)) # will take the whole screen
x = 1:70; y = 1:70 # make variables 
x.val = (x^2);   y.val = y
xy.out = outer(x.val, y.val, "-") # manipulate the variables
nb.col.pal = length(grDevices::hcl.pals()) # find all palettes in hcl.pals
rdm.nb = sample(c(1:nb.col.pal), size = 1); print(rdm.nb) # Get random number
# [1] 59
xy.out = xy.out + rnorm(length(xy.out), mean = 0, sd = 200) # Prettify with RDM
image(t(xy.out),xaxt = "n", yaxt = "n", bty = "n",
      col = hcl.colors(1000, palette = hcl.pals()[rdm.nb])) # 61,93 looks good

Challenge 1

Select randomly 4 values out of a vector of numbers ranging from 1 to 10 : 1. without replacement and 2. with replacement.

Use sample() to perform this task.

sample()

Challenge 1 - Solution

Select randomly 4 values out of a vector of numbers ranging from 1 to 10 : 1. without replacement and 2. with replacement.

set.seed(12345) # Sets the random number generator to a fix value
vec.1.10 = 1:10 # Make the vector to choose from 
sample(x = vec.1.10, size = 4, replace = FALSE) # Sample 4 nb without replacement
# [1] 3 8 2 5
sample(x = vec.1.10, size = 4, replace = TRUE) # Sample 4 nb with replacement
# [1] 8 2 6 6

As you can see in the last example, there are 2 "6"s, since each time a random number was picked, all numbers could be randomly chosen.

1000 draws with replacement from 0 or 1, with equal probability for each.

set.seed(123); table(sample(x = 0:1, size = 1000, replace = T,prob = c(.5,.5)))
# 
#   0   1 
# 493 507

Challenge 2

Create a data frame with variable

• x ranging from 1 to 10,
• $y = 2+ 3 * x$ and
• a grouping factor (gr) with one group with the 5 smallest values and another group with the 5 largest values.

Challenge 2 - Solution

x = 1:10
y = 2 + 3 * x
gr = rep(letters[1:2],each = 5)
linear.df = data.frame(x,y,gr)
plot(y~x, col = as.factor(gr), data = linear.df, pch = 19, cex = 1.5)

Are you sure?

• What characterizes a random phenomenon or all situations that involves chance or a likelihood that certain outcomes are realized?
• A degree of uncertainty!
• Below are 3 and 100 coins (the lines) flipped 10 and 5000 times each (the x-axis).

Probability and odds

• When describing the "chance" or "likelihood" of an outcome, you can use a probability or odds (which reports probabilities in terms of ratios)
• They are related:
  • $p_{success} = \frac{\text{Nb of successes}}{\text{Total number of events}}$ and $1-p_{success} \text{ or } q = \frac{\text{Nb of failures}}{\text{Total number of events}}$
  • odds in favor of an event occurring is $odds_{favour} = p/(1-p)$ or $odds = \frac{\text{probability of success}}{\text{probability of failure}}$ or $odds_{success} = \frac{\frac{\text{Nb of successes}}{\text{Total number of events}}}{\frac{\text{Nb of failures}}{\text{Total number of events}}}=\frac{\text{Nb of successes}}{\text{Nb of failures}}$. Otherwise, the odds against is $odds_{against} = (1-p)/p$
• Why bother using probability or odds?
• probability ranges between 0 and 1, but odds range from 0 to $\infty$.

• It is useful to change the scale characterizing the chance of an event occurring. But from the "odds scale" on the right, the odds of losing is not symmetric (between 0 and 1) compared to odds of winning (between 1 and $\infty$). How can we solve that and make both sides symmetric?

Probability and odds

• By using the log of the odds.
• Log (natural log, or ln) odds can now
  • range from $-\infty$ to $\infty$ and
  • makes the (log)odds-scale symmetric.
  • It is calculated as $ln (odds) = ln(p/(1-p))$ (i.e., the logit function).
• This will become handy when building a logistic model as this is the 'logit' link

Probability and odds

• Example: a horse runs 100 races and wins 20 of them. What is the odds of winning?
• $odds = p/(1-p)$ or $odds = \frac{\text{probability of success}}{\text{probability of failure}}$

tot.nb.ev = 100; success = 20
failure = tot.nb.ev-success
p = success/tot.nb.ev; q = failure/tot.nb.ev
(odds.favor = (success)/(failure)) # Odds in favor of event ( 1 in 4), same as p/q
# [1] 0.25
(odds.agnst = (failure)/(success)) # Odds against the event (here 4 to 1), same as q/p
# [1] 4

The horse wins 1 race to 4 fails (odds of winning = 0.25). So, for 5 races it will win 1 and loose 4.

Conversely, the horse odds of failing is 4:1 (read 4 to 1). The horse is 4 times more likely to fail than to succeed in a race. Odds can also be represented as ratios (and integer ratios).

Probability and odds

• Example: a horse runs 100 races and wins 20 of them. What is the probability of winning?
• $p = odds/(1+odds)$

tot.nb.ev = 100; success = 20
failure = tot.nb.ev-success
(probability.favor = (success)/(tot.nb.ev)) # probability of event (here 20%), fractions(probability.favor)
# [1] 0.2
(q = failure/tot.nb.ev)
# [1] 0.8

• So the horse wins 1/5 or 20% of the times.
• The probability of a particular outcome (say getting a head for a coin) for a particular event (tossing a coin) is the number of a particular outcome divided by the total number of results from the event.
• You can see that $odds = p/(1-p) = (.2)/(1-.2) = .2/.8 = 1/4$

Probability and odds

• Let's code the difference between probability, odds and log odds.

# Get some LOG ODDS numbers 
log_odds = seq(from = -5, to = 5, by = .25)
odds = exp(log_odds) # Transformed into odds 
# Make inverse logit function 
inv.logit = function(x) {exp(x)/(1 + exp(x))} # takes log odds as input. Same as function(x){1/(1 + exp(-x))}
p = inv.logit(log_odds) # This is  p = odds/(1 + odds) 
q = 1-p # Probability of failure (1-p)
# Store log_odds other data to plot 
d = data.frame(log_odds, odds, p, q) 
format(x = ht(d,3),digits = 2, scientific = F)
#         log_odds     odds      p      q
# head.1      -5.0   0.0067 0.0067 0.9933
# head.2      -4.8   0.0087 0.0086 0.9914
# head.3      -4.5   0.0111 0.0110 0.9890
# tail.39      4.5  90.0171 0.9890 0.0110
# tail.40      4.8 115.5843 0.9914 0.0086
# tail.41      5.0 148.4132 0.9933 0.0067

Probability and odds

• The three plots show the contrast between probability, odds and log odds.
• The regions in blue (success) and red (failure) are highlighting the part of the function that is considering a success and a failure for the black line.

Probability and p-values

• What is the difference between a probability and a p-value?
• A p-value is not only the probability of a specific event.
• A p-value consists of 3 parts. They are the :
  1. Probability of an event 'A' or that random chance generated the data 'A' (or would give the result we got) PLUS
  2. Probability of observing something else (other events) that is equally rare PLUS
  3. Probability of observing something else that is rarer (or have less chance to occur or more extreme) than the event 'A'.
• NOTE 1: a p-value for continuous variables is the area under the curve AT the value and GREATER or LOWER the CRITICAL value or BOTH (greater and lower) depending on the test.
• NOTE 2: the probability of an event occurring is not the same as the p-value (sum of the 3 parts).
  • We add part 2. and 3. because it gives us information about the specialness of this rare event.
  • In other words, we want to know if the event we measured is statistically special or, as you might have heard it, statistically significant!

Probability and p-values

What is the p-value 1. for getting a human height of 108.04 cm? 2. What is the probability of having a height between 109.9995 and 110.0005 cm? What is the p-value for someone of height 110 cm?. Below are the distributions of hypothetical heights with mean = 110 and sd = 1.

#         lb       ub
# 1 109.9995 110.0005

The answer is 1. 0.05 or 5%, 2. 0.04%, 3. 1 or 100%!

Probability distributions are crucial for simulations

• When performing simulations, one has to keep in mind the
  1. type of distribution underlying the data of interest
  2. "arguments" or "parameters" of the distribution itself (mean, standard deviation, rate, degrees of freedom, etc.)
  3. statistical model making the relationship between the response variable and the explanatory variable.
  4. ideas or things we want to show or learn with the simulation

Why statistical distributions matter?

• As we discussed in the introduction, we are interested in what is underlying natural processes.
• To mimic the natural processes,
  • we need a way to get some similar data as we would if we collect them (we need a way to generate (pseudo-)random variables).
  • In other words, to help us in performing simulations, we need to get some data (random variables, X) which numerical values are outcomes of random processes.
• The idea here is to find the distribution underlying a random process which generate a particular phenomenon in nature and use it to output values to be used in your simulation.
• After generating the data, you can analyze it and see how it can be manipulated.

What is a random variable?

• Imagine a number is presented to you like so:

# [1] 13.9

• From which distribution is this variable coming from?
• Would be helpful to have more information

What is a random variable?

• Now imagine you are able to get more numbers from the same process:

# [1] 13.9 14.5 18.1 15.1 15.3

• This is more informative. We can take the mean and the standard deviation of these numbers: mean = 15.4, sd = 1.6.
• Again, from which distribution is this variable coming from?
• We could plot the result to get a feel of the distribution:

• Still hard to definitely tell

What is a random variable?

• What if you get even more numbers (say 100) from the same random process:

#   [1] 13.9 14.5 18.1 15.1 15.3 18.4 15.9 12.5 13.6 14.1 17.4 15.7 15.8 15.2 13.9
#  [16] 18.6 16.0 11.1 16.4 14.1 12.9 14.6 12.9 13.5 13.7 11.6 16.7 15.3 12.7 17.5
#  [31] 15.9 14.4 16.8 16.8 16.6 16.4 16.1 14.9 14.4 14.2 13.6 14.6 12.5 19.3 17.4
#  [46] 12.8 14.2 14.1 16.6 14.8 15.5 14.9 14.9 17.7 14.5 18.0 11.9 16.2 15.2 15.4
#  [61] 15.8 14.0 14.3 13.0 12.9 15.6 15.9 15.1 16.8 19.1 14.0 10.4 17.0 13.6 13.6
#  [76] 17.1 14.4 12.6 15.4 14.7 15.0 15.8 14.3 16.3 14.6 15.7 17.2 15.9 14.3 17.3
#  [91] 17.0 16.1 15.5 13.7 17.7 13.8 19.4 18.1 14.5 12.9

• The mean = 15.4, sd = 1.6.
• Check the graph

• That looks normal!

Functions for various distributions

• In R, there is a convention on how to call different distributions.

• In the next slides, you'll see some built in functions to use these distributions

Functions for various distributions

• Here is an example of the d, p, r, q functions for the normal distribution.
• Note that lower.tail = TRUE by default. It means that the probability is calculated from the "left" whereas lower.tail = FALSE means that the probability will be calculated from the "right".

Functions for various distributions (special)

# probability of > 2 people with the same birthday
(a = pbirthday(23, coincident = 3))
# [1] 0.01441541
# 0.9 probability of >=3 coincident birthdays
(a = qbirthday(coincident = 3, prob = 0.9)) # Gives the number of people
# [1] 135
# Chance of >=4  coincident birthdays in 300 people
(a = pbirthday(300, coincident = 4))
# [1] 0.9740397

Functions for various distributions

Functions for various distributions

Common statistical distributions (density or mass)

Summary of important distribution

• Distributions of random variables comes in 2 flavors:
  1. Discrete (infinite integer number of distinct values, e.g., 0, 5, 19,...)
  2. Continuous (infinite decimal number of possible values, e.g., 0.1,0.62113,...)

• see Probability Cheat Sheet from P. Michaeli, from J. Fieberg, or from W. Chen for more details

Selecting the right distribution

• See "Figure 6A.15 from Probabilistic approaches to risk by Aswath Damodaran"

The same distribution?

Careful!

• Some distribution are similar but don't show the same process

Almost all roads lead to "the normal distribution"

The image below shows a small roadmap demonstrating the link between different statistical distributions.

Other distributions with packages

• Other packages have other distributions. See the CRAN Task View: Probability Distributions if you need a funky distribution that doesn't come in base R.
• Some distributions can be recreated in R even if they haven't their own function name (e.g., Bernoulli).
  • Bernoulli distribution can be simulated with a binomial distribution with size = 1.

Binomial distribution

• Binomial distribution is finite (number of trials is not infinite) and discrete.

• It assumes that:

  • Each independent trial is either a success (1, T, Yes) or failure (0, F, No)
  • There is a fixed number of trials and the probability on each trial is constant.

• The probability mass function is:

$$f(x) = {n \choose x} p^x(1-p)^{n-x}$$ where

• $p$ is the probability of success in a single trial ($(1-p)=q$ is the probability of failure)
• $n$ is the number of trials
• $x$ is the number of trials that are successes
• ${n \choose x} = \frac{n!}{r!(n-r)!}$ is the number of combinations where the order doesn't matter.

Binomial distribution

• Examples
  • Number of heads in a row after 3 tosses (heads or tail)
  • Wright-Fisher model (Genetic drift, p and q)
  • Testing the side effects of a new medications (with and without effect)
  • Model river overflows due to excessive rain (overflows or not).
  • Number of individuals that have a survival probability in an experiment.

Binomial distribution in R

• In R, you can generate binomial distributed variables with the rbinom function.

rbinom(n = number of repetitions, # or number of 'experiments' or observations
       size = sample size, # between 0 and up to the size (0:size)
       p = probability of success)

• Think of it as the number of times you'd get a success, with probability $p$, for a number of times you repeat the experiment.

set.seed(12345)
# 1000 experiments where each time, e.g., flipping a coin, I can either have a
table(rbinom(n = 1000, size=1, prob=0.5), dnn=NULL) # success or failure with p=.5
#   0   1 
# 469 531
# 1 experiment with 1000 coins summing all successes with p=.5
rbinom(n = 1, size=1000, prob=0.5)
# [1] 476
# 1000 experiments where each time, for example flipping 10 coins, where I 
rbind(random = table(rbinom(n = 1000, size=10, prob=0.5), dnn=NULL)/1000,# sum the success with p=.5
      theory = round(dbinom(x = 0:10, size=10, prob=0.5),3))
#            0     1     2     3     4     5     6     7     8     9    10
# random 0.001 0.009 0.045 0.105 0.191 0.260 0.211 0.107 0.058 0.012 0.001
# theory 0.001 0.010 0.044 0.117 0.205 0.246 0.205 0.117 0.044 0.010 0.001

Binomial distribution

• Example 1: probability of having a certain number of offspring with a certain genotype (e.g., Aa) given the genotypes of the parents
  • If a male (AA) mates a female (Aa), what is the probability that 6 of their 7 offspring is Aa (7 choose 6)
  • The probability of getting Aa from AA x Aa is [1/2 AA and 1/2 Aa],
  • $P(6|7) = {7 \choose 6} (1/2)^6(1-1/2)^{7-6} = 7/128$ (about 0.055)

dbinom(x = 6,7,prob = .5)
# [1] 0.0546875

• Example 2: Probability of having 3 daughters if you have 5 children.
  • if you have 5 children, these are independent trials, so n = 5
  • "success" is the probability that the child is a daughter (p = 1/2)
  • $P(3|5) = {5 \choose 3} (1/2)^3(1-1/2)^{5-3} = 5/16$ (about 0.31)

dbinom(x = 3,5,prob = .5) 
# [1] 0.3125

Binomial distribution

• This is a visualization of the previous examples.

plot(0:8, dbinom(x = 0:8, size=7, prob=0.5),type='h', col = col.b, lwd=lwd, ylim = c(0,.4), ylab = "Probability", main = "s=7 (blue) and s=5 (red); p=0.5 "); add.l(.05)
points(x = 6, y = dbinom(x = 6, size=7, prob=0.5), pch =19, col = col.b)
lines(c(0:8)+0.1, dbinom(x = 0:8, size=5, prob=0.5),type='h', col = col.r, lwd=lwd)
points(x = 3+.1, y = dbinom(x = 3, size=5, prob=0.5), pch =19, col = col.r)
legend("topright",legend = c("Genetics","Daughter"), lty =c(1), col = c("blue","red"), lwd = 3)

Binomial distribution

# For dbinom, x is the vector of quantiles 
plot(0:1, dbinom(x = 0:1, size=1, prob=0.5),type='h', col=col.b, lwd=lwd, ylim=c(0,1), ylab="Probability", main = "s=1, p=.5"); add.l()
plot(0:20, dbinom(x=0:20, size=20, prob=0.4), type='h', col=col.b, lwd=lwd, ylim=c(0,.2), ylab="", main = "s=20, p=.4");add.l()
plot(0:23, dbinom(x=0:23, size=20, prob=0.9), type='h', col=col.b, lwd=lwd, ylim=c(0,0.3), ylab="", main = "s=20, p=.9");add.l()
plot(40:110, dbinom(x=40:110, size=150, prob=0.5), type='h', col=col.b, lwd=lwd, ylim=c(0,.15), ylab="Probability", main="s=150, p=.5");add.l()
plot(40:51, dbinom(x = 40:51, size=50, prob=0.99),type='h', col = col.b, lwd=lwd, ylim = c(0,.8), ylab = "", main = "s=50, p=.99");add.l()
hist(rbinom(n = 100000, size=50, prob=0.99),xlab = "Nb of successes", breaks = 100, col = col.b, main = "s=50, p=.99", probability = F)

Challenge 3

What is the probability (in %), of getting 4 heads and 4 tails from 4 tosses of a fair coin? Use the binomial distribution to answer this question.

Challenge 3 - Solution

What is the probability (in %), of getting 4 heads or 4 tails from 4 tosses of a fair coin?

dbinom(x = 0, size = 4, prob = .5) + dbinom(x = 4, size = 4, prob = .5)
# [1] 0.125
curve(dbinom(x,4,.5),0,6, n = 7,type = "h", lwd = lwd)
curve(dbinom(x,4,.5),0,0, type = "h",add = TRUE, col = "red", lwd = lwd)
curve(dbinom(x,4,.5),4,4, type = "h",add = TRUE, col = "red", lwd = lwd)

Poisson distribution

• To be Poisson distributed, the data generally assumes that

  • the probability of some event is small over a short period of time (or area), but that there are a lot of events (as there would be a lot of time passing, or a lot of space) in which the process could happen
  • The number of event reported is independent from one time point to the other.
  • In summary, it models rare events that happen at a certain rate. In other words, there are many opportunities to succeed (n is large), but the probability of success (for each trial) is small.

• The probability mass function is:

$$f(x) = e^{-\lambda}\frac{\lambda^x}{x!}$$ where $\mu = \sigma^2 = \lambda$

Poisson distribution

• Examples (something per unit time or space for example):
  • Nb of deaths due to a disease (e.g., SARS-CoV 2) for a period of time
  • Ticks of a radiation counter (Geiger counter which detects and measures ionizing radiation)
  • DNA mutation
  • Colonization of animals on remote islands
  • Mutation models
  • Recombination models
  • Bacteria colonies growing on an agar plate

Poisson distribution in R

• In R, you can generate Poisson distributed variables with the rpois function. Below is the structure of the function.

rpois(n = number of repetitions,
      lambda = rate) # doesn't need to be an integer

• You can think of this function as the number of times you'd get a success, when the number of events per unit of time or space is lambda
• E.g., with a camera trap, I see on average 2 wolves/day (lambda is this rate), then if I go on the field for one day, what could be the number of wolves that I observe?

set.seed(5937); rpois(n = 1,lambda = 2)
# [1] 5

• Wow! That day was amazing! The probability of seeing this number is 0.04 and the probability of seeing 5 or more is 0.02

Poisson distribution

add.l <- function(by=.1) {  abline(h=seq(0,1, by = by), lty = 3, lwd = .3)}
plot(0:5, dpois(x = 0:5, lambda=0.05),type='h', col = col.b, lwd=lwd, ylim = c(0,1), ylab = "Probability", main = "lambda = 0.05");add.l();abline(v = 0.05, lty = 3)
plot(0:10, dpois(x = 0:10, lambda=1),type='h', col = col.b, lwd=lwd, ylim = c(0,1), ylab = "Probability", main = "lambda = 1"); add.l();abline(v = 1, lty = 3)
plot(0:20, dpois(x = 0:20, lambda=5),type='h', col = col.b, lwd=lwd, ylim = c(0,0.3), ylab = "Probability", main = "lambda = 5");add.l();abline(v = 5, lty = 3)
plot(0:40, dpois(x = 0:40, lambda=20),type='h', col = col.b, lwd=lwd, ylim = c(0,.2), ylab = "Probability", main = "lambda = 20");add.l();abline(v = 20, lty = 3)
hist(rpois(n = 1000, lambda=25), xlab = "Nb of successes", breaks = 100, col = col.b, main = "n = 1000, lambda = 25", probability = F);abline(v = 25, lty = 3)
hist(rpois(n = 100000, lambda=50),xlab = "Nb of successes", breaks = 100, col = col.b, main = "n = 100000, lambda = 50", probability = F);abline(v = 50, lty = 3)

Poisson distribution vs binomial

• Difference between Binomial and Poisson distributions:
  • Number of events in Poisson can be infinite (so n tends towards infinity and the probability of each event approaches 0)
  • Number of trials is finite ( $n$ ) in Binomial distribution

# [1] 0.6766764
# [1] 0.6766764

Uniform distribution (discrete)

• This distribution can take any discrete value in a range

# Define uniform discrete 
dunifdisc<-function(x, min=0, max=1) ifelse(x>=min & x<=max & round(x)==x, 1/(max-min+1), 0)
runifdisc<-function(n, min=0, max=1) sample(min:max, n, replace=T)
curve(dunifdisc(x, 7,10), type = "h",from=6, to=11, col="black", xlab = "",
      lwd=1, ylim = c(0,.5), ylab = "Probability", main = "Uniform discrete"); title(xlab = "x",line=2.2, cex.lab=1.2)

set.seed(12345)
sample(7:10, size = 1,replace = T)
# [1] 8
runifdisc(1,7,10)
# [1] 9

Uniform distribution (continuous)

• This distribution will give any (continuous) value between a minimum and a maximum value.
• This can come handy when generating angles.

x = seq(0,360,by = .1)
y <- dunif(x, min = min(x),max = max(x))
plot(y~x, type = "n", ylim = c(0,1.5*max(y)), ylab = "Density", xlab = "", main = "Uniform continuous")
polygon(c(x, rev(x), 0), c(y, rep(0,length(y)), 0), col=scales::alpha("blue",.5))
title(xlab = "x",line=2.2, cex.lab=1.2)

set.seed(12345)
runif(n = 1, min = 0, max = 360)
# [1] 259.5254

Normal distribution (Gaussian)

• This distribution is the most important to know (see Central limit theorem).
• It has interesting properties (mean: center, sd: spread, bell-shaped, symmetry, etc.)
• Can approximate other distribution (might require data transformation)
• The normal distribution with mean = 0 and sd = 1 is called the Standard Normal

The general form of the probability density function for a normal distribution is

$$f(x)={\frac{1}{\sigma\sqrt{2\pi}}}e^{-{\frac{1}{2}}\Bigl(\frac {x-\mu}{\sigma}\Bigr)^2}$$

The important part here is the mean $\mu$ and the standard deviation $\sigma$. Usually summarized as $N(\mu, \sigma)$. However, note that sometimes, $\sigma$ is the variance, but in R, it is the standard deviation.

Here is a link to play with a Normal distribution interactively.

Normal distribution (Gaussian)

• Examples:
  • Distribution of a lot of continuous phenotypes
  • Arm length, or human height in a population
  • velocity of a collection of molecules in a gas of liquid
  • the error from measurements (length of beak).

Normal distribution (Gaussian)

• The iris dataset contains values that are 'drawn' from a normal distribution.
• Here are the histogram for sepal length (mm).

Normal distribution (dissection)

curve(expr = dnorm(x = x, mean=0,sd=1), 
      from = -5, to = 5, ylim = c(0,1),
      col="black", ylab = "Density",
      main = "Density normal") 
abline(h=seq(0,1, by = .1), 
       lty = 3, lwd = .3)

Convert a Normal dist. to Standard Normal

• If your data is normally distributed (i.e., $X \sim N(\mu, \sigma)$), you can convert the data to be $X \sim N(0, 1)$ with a scaling $$Z = \frac{(X-\mu)}{\sigma}$$

Quantile value for a normal distribution

If X is a random variable distributed with a Standard normal distribution, what is the probability of finding X less or equal to $1.645$. In other words, $P(X\leq 1.645)$ where $X\sim N(0,1)$?

# Gives the probability at X 
pnorm(1.645) 
# [1] 0.9500151
# Gives X at a prob.
qnorm(p = 0.05, lower.tail = F) 
# [1] 1.644854
# Value for 5% in both tails
qnorm(p = 0.025, lower.tail = F) 
# [1] 1.959964

Normal distribution

# The function is not shown, but can be found in the script (Markdown)
draw.normal(where = "both",  prob = 0.05/2) # 2.5% in each tail
draw.normal(where = "both",  prob = 0.2/2 ) # 10% in each tail 
draw.normal(where = "both",  prob = 0.5/2 ) # 25% in each tail 
draw.normal(where = "both",  prob = 0.95/2) # 47.5 % in each tail
draw.normal(where = "left",  prob = 0.05  ) # 5% only in the left tail 
draw.normal(where = "right", prob = 0.05  ) # 5% only in the right tail

Normal distribution (pnorm)

Get the probability of finding data relative to a standard deviation number

sd = 1
probability.left.side  = (pnorm(q = c(sd*1, sd*2, sd*3), lower.tail = F))
probability.right.side = (pnorm(q = c(sd*1, sd*2, sd*3), lower.tail = T))
percent.data.under.curve = probability.right.side - probability.left.side
p.from.mean = round(percent.data.under.curve*100,2)

• So from the mean of the standard normal distribution:
  • if you are at $\mu \pm 1 \sigma$ , you have 68.27% of the data.
  • At $\mu \pm 2 \sigma$ it's 95.45% and at $\mu \pm 3 \sigma$ it's 99.73%

Normal distribution (qnorm)

Find the "x" value (quantile) based on the probability (area under the curve between 2 values) of a standard normal distribution.

qnorm(p = c(.75, .95,.975, .995), mean = 0, sd = 1, lower.tail = F) # notice the "lower.tail"
# [1] -0.6744898 -1.6448536 -1.9599640 -2.5758293
qnorm(p = c(.75, .95,.975, .995), mean = 0, sd = 1, lower.tail = T)
# [1] 0.6744898 1.6448536 1.9599640 2.5758293

• for a probability of $50\%$ between the "x" value is -0.67, 0.67.
• for a probability of $90\%$ between the "x" value is -1.64, 1.64.
• for a probability of $95\%$ between the "x" value is -1.96, 1.96.
• for a probability of $99\%$ between the "x" value is -2.58, 2.58.

Normal distribution

• Hypothesis testing: when we are trying to find differences in the mean of a sample and a theoretical value, this is what we are assuming.

Normal distribution

• Let's see what happens if we increase the sample size from 10 to 100.

Challenge 4

Let's say you want to simulate the length of different beak sizes for 10 birds.

Use rnorm() to generate 10 random normal numbers.

rnorm()

Challenge 4 - Solution

Let's say you want to simulate the length of different beak sizes for 10 birds.

Use rnorm() to generate 10 random normal numbers.

set.seed(1234)
n <-10
rnorm(n)
#  [1] -1.2070657  0.2774292  1.0844412 -2.3456977  0.4291247  0.5060559
#  [7] -0.5747400 -0.5466319 -0.5644520 -0.8900378

• Why are there negative values? Because the mean = 0, so some values spill over the mean on both sides, and in the negative numbers.

Challenge 5

What is the probability of colonizing an island if a species has a 1 in a million (1/1000000) chance of colonizing it and that we wait 1000000 (1 million) years?

Hint: think of colonization as a success.

• Think of the problem for ~1 min;
• Pair with a neighbor and hare your ideas and code a solution (~4 min)
• Discuss with the group your findings (~5 min).

Challenge 5 - Solution

size = 1000000
prob = 1/1000000
x = seq(0,10,by = 1)
plot(dbinom(x = x, size = size, prob = prob)~x, 
     type = "h", lwd = 4, ylab = "Probability", xlab = "Event")
points(dbinom(x = x[-1], size = size, prob = prob)~x[-1], 
       type = "h", lwd = 4, col = "red")

# This is the probability 
pbinom(0,size = size,prob = prob,lower.tail = F)
# [1] 0.6321207
# pbinom(0,size = size,prob = prob,lower.tail = T)

Central limit theorem (CLT)

• It has been said that if "physicists have general relativity, chemists have the periodic table of elements and biologists have evolution or DNA, statisticians have the central limit theorem" (Nathan Uyttendaele).

• The Central limit theorem (CLT) stipulates that:

  • The more you add (sum) identically distributed random variables, the more the distribution function of the new variable will converge to a normal distribution (or could be approximated by a normal distribution).

• We will use simulations to test this assertion!

Challenge CLT

• Imagine I give you a set of dice to conduct an experiment in which you are going to roll them 10,000 times.
• X is a random variable which represents the sum of the 2 numbers after rolling the die
  • $S = {(x_1, x_2, ..., x_n): x_1 = 1,2,3,4,5,6; x_2 = 1,2,3,4,5,6; ...; x_n = 1,2,3,4,5,6}$
  • $X(\omega) = x_1+x_2+...+x_n$ for $\omega = (x_1, x_2, ..., x_n) \in S$

Make a histogram of the sum of the values from rolling a set of 1, 2, 10 and 50 die 10,000 times.

Hint: to roll a dice, use the following code

set.seed(42)
sample(1:6, size = 1, replace = TRUE)

• Note: these are discrete values

Challenge CLT - Solution

Make a histogram of the sum of the values from rolling a set of 1, 2, 10 and 50 die 10,000 times.

set.seed(12345)
roll = 1e4 # Nb of times we do the experiment (characterize the underlying r.v.)
for (nb.dice in c(1,2,10,50)) {
  res.exp = replicate(roll, simplify = T,
                      sample(1:6, size = nb.dice, replace = TRUE))
  if(nb.dice == 1){Xsum = res.exp } 
  else {Xsum = apply(res.exp, 2, sum)} # add independent r.v. (sum of cols)
  hist(Xsum, main=paste("n =",roll,ifelse(nb.dice==1,"dice =","die ="),nb.dice))}

Central limit theorem (CLT)

• The CLT could be shown by adding more and more random variables distributed as an exponential with rate = 1.

n = 1000 # Number of points 
# Generate multiple additions of random variables 
for(i in c(2, 50, 1000, 5000)){
  clt = replicate(i, rexp(n, rate = 1), simplify = FALSE)
  hist(apply(do.call(cbind,clt),1,sum), main = paste("Hist. of",i,"variables"), xlab = "x") # Draw the histogram 
}

Central limit theorem (CLT)

• Same thing as above, but for a random variable coming from a uniform discrete distribution (with min = 0, max = 1)!

# Generate multiple additions of random variables 
for(i in c(2, 50, 1000, 5000)){
  clt = replicate(i, runifdisc(n),simplify = FALSE)
  hist(apply(do.call(cbind,clt),1,sum), main = paste("Histogram of",i,"variables"),xlab = "x") # Draw the histogram 
}

Central limit theorem (CLT)

• You are asked to throw 20 dice. If the sum of all 20 dice is $\geq 100$, you get a candy, if not, you don't get one. What is the probability of having the candy?

n = 1e6 # Number of points 
# Generate multiple additions of random variables 
clt = replicate(n = 20, # flip 20 dice
                runifdisc(n = n, min = 1, max = 6), simplify = FALSE)
sum.rdm.var = apply(do.call(what = cbind, args = clt), 1, sum) 
mean.rdmv = mean(sum.rdm.var); sd.rdmv = sd(sum.rdm.var) # mean and sd
hist(sum.rdm.var, main = paste("Histogram of",20,"variables"), xlab = "x")
curve(expr = dnorm(x,mean = mean.rdmv,sd = sd.rdmv),
      from = mean.rdmv-5*sd.rdmv, to = mean.rdmv+5*sd.rdmv,ylab = "Density")

• Using the CLT, the probability of getting a sum $\geq 100$ is approx. 4.36e-05.

Central limit theorem (CLT)

• Galton's board (quincunx) is a representation of a lot of Bernoulli trials that will form a normal distribution

• Below is a simulation of 1000 balls, hitting 500 pegs that were left from the point "0". The normal curve was drawn with a standard deviation of $\sqrt{n}-1$

Multivariate normal distribution

• There is a way to generate multiple normal distributions at once from a covariance matrix.

set.seed(24601) # setting this so the random results will be repeatable 
library(MASS)
covmat <- matrix(c(1.0,   0.2,   0.6, # variance covariance matrix of the data 
                   0.2,   2.0,  -0.5, # Positive-definite symmetric matrix
                   0.6,  -0.5,   1.0), nrow=3) 
data <- mvrnorm(n = 300,  
                mu = c(1,-1,0), # mean of the data 
                Sigma=covmat) # generate random data that match that variance covariance matrix
plot(data[,1:2], pch = 19, col = b.5); abline(h=0,v=0,lty = 3)
plot(data[,2:3], pch = 19, col = b.5); abline(h=0,v=0,lty = 3)

Multivariate normal distribution

• There is a way to generate multiple normal distributions at once from a covariance matrix.

ggplot2::autoplot(vegan::rda(data)) + theme_classic() + 
  theme(plot.margin = unit(c(0,0.5,0,0.5), "cm"), legend.position="none")+
  geom_hline(yintercept=0, linetype="dotted",color = "grey50", size=0.5)+
  geom_vline(xintercept=0, linetype="dotted",color = "grey50", size=0.5)

Multivariate normal distribution

• You can also generate multiple normal distributions at once with a correlation matrix between the variables to get the covariance matrix.
• The diagonal are variances for each variable, the upper and lower triangles are the same and refer to the covariance between the variables. The faux package can help you simulate factorial designs and multiple normal correlated variables (see rnorm_multi).

library(faux); set.seed(12345); n = 500 #number of values to generate for the variable x
traits = rnorm_multi(n = n, # ??rnorm_multi
                     mu = c(11.317531, 10.470744, 9.377967), # Mean of 3 variables
                     sd = c( 0.6247866, 0.7134755, 0.5502800), # SD of variables 
                     r = c(0.4887644, 0.4678446, 0.7056161), # correlation between variables 
                     varnames = c("Beak length", "Beak depth", "Beak width"), empirical = TRUE)
plot(traits[,1:2],col=b.5,pch=19,main="Simulated data",xlab="");title(xlab="Beak length",ylab="",line=2.2,cex.lab=1.2)

Multivariate normal: PCA simulation

set.seed(123) # setting this so the random results will be repeatable 
library(MASS)
# Simulating 3 traits for 4 different species 
n = 200
Amat1 = MASS::mvrnorm(n, mu = c(11.2,11.8,9.91), Sigma = diag(c(1.31,1.01,1.02))) 
Amat2 = MASS::mvrnorm(n, mu = c(7.16,8.54,6.82), Sigma = diag(c(0.445,0.546,0.350)))
Amat3 = MASS::mvrnorm(n, mu = c(15.6,14.6,13.5), Sigma = diag(c(1.43,0.885,0.990)))
Amat4 = MASS::mvrnorm(n, mu = c(8.65,14.1,8.24), Sigma = diag(c(0.535,0.844,0.426)))
Amat = rbind(Amat1,Amat2,Amat3,Amat4)
Amat.gr = cbind(Amat, gl(4,k=n,labels = c(1,2,3,4)))
# by(Amat.gr[,1:3],INDICES = Amat.gr[,4],FUN = cov) # calc. cov. mat. for all gr 
summary(m1 <- prcomp(Amat, scale= T))
# Importance of components:
#                           PC1    PC2    PC3
# Standard deviation     1.5566 0.6977 0.3005
# Proportion of Variance 0.8076 0.1623 0.0301
# Cumulative Proportion  0.8076 0.9699 1.0000

Multivariate normal: PCA simulation

# biplot(m1, xlabs=rep(".", nrow(Amat)), cex = 3)
plot(Amat[,1],Amat[,2], pch = 19, col = gl(4,k=n,labels = c(1,2,3,4))); abline(h=mean(Amat[,1]),v=mean(Amat[,2]),lty = 3)
plot(vegan::scores(m1), asp = 1, pch = 19, col = gl(4,k=n,labels = c(1,2,3,4))); abline(h=0,v=0,lty = 3)

# library(ggvegan)
# autoplot(vegan::rda(Amat))

Find the parameters of a distribution

x = rnorm(100)
# x = rt(100,1000)
library(fitdistrplus)
descdist(x, discrete = FALSE)

# summary statistics
# ------
# min:  -2.307474   max:  2.735209 
# median:  -0.06346588 
# mean:  -0.03805934 
# estimated sd:  0.9945627 
# estimated skewness:  0.2436371 
# estimated kurtosis:  2.934408
fit.pois <- fitdist(x, "norm")
plot(fit.pois)

Linear model (lm) refresher

$$Y = \beta_{0} + \beta_{1} x_{1} + \cdots + \beta_{p} x_{p} + \epsilon$$

• $Y$ is the response variable
• $\beta_0$ is the intercept
• $\beta_1$ is the coefficient of variation for the first explanatory variable ( $x_1$ )
• $\beta_p$ is the coefficient of variation for the $p^{th}$ explanatory variable for the $p^{th}$ $x_p$ explanatory variable
• $\epsilon$ is the residual of the model. Note that $\epsilon \sim N(\mu=0,sd = \sigma)$

Linear model (lm) refresher

• The goal is to find the best estimation of the parameters ( $\beta$ s), while minimizing the residuals, and assess the goodness of fit of the model
• We use the least-squares in order to fit the line to the data.
• In order to interpret your results correctly when using a linear model, be sure you are not crossing the LINE:
  • L: Linearity, — I: Independence, — N: normally distributed residuals, — E: Equality of variance

Linear model (lm) refresher

• You can add:
  1. the confidence limits for the expected values (confidence interval or uncertainty around the mean predictions (y) for a value of x given the parameters of the model; blue shaded polygon predict(lm.out, interval = "confidence"); at mean of x and y, the uncertainty comes only from the intercept) to the line and
  2. the 95% central quantiles of the Normal distribution (or prediction interval, grey shaded polygon for the simulated one; compared with red shading lines which is the theoretical one predict(lm.out, interval="prediction"))

Generalized linear model (GLM) refresher

• What we've seen so far are cases in which the response variable was continuous and where the "error" is normally distributed.
• But there might be cases where the response variable is not continuous or has an error structure that is not normal.
• Therefore, we need to look for GLMs.

• Keep in mind that the x-values don't need to be normally distributed.

Generalized linear model (GLM) refresher

• As we have seen, the idea being linear regression is to model the expected value. We are doing the same thing for a poisson regression model.

Simulate categories (t-test, Anova)

set.seed(1234); n = 1000
y1 = rnorm(n, mean = 15, sd = 1)
y2 = rnorm(n, mean = 15.5, sd = 1)
sim.aov1 = data.frame(y = y1, gr = "A")
sim.aov2 = data.frame(y = y2, gr = "B")
df.aov = rbind(sim.aov1, sim.aov2)
df.aov$gr = factor(df.aov$gr)
# t.test(y~gr, data = df.aov) or
aov.out = aov(y~gr, data = df.aov)
#summary(aov.out)
tk.test = TukeyHSD(aov.out)
round(tk.test$gr,2)
#     diff  lwr  upr p adj
# B-A 0.54 0.45 0.63     0

plot(y~gr, data = df.aov)

GLM: Linear model

• We start by generating a linear model with an intercept and slope for 100 points (here random uniform between 0 and 1), an the error in the measurements.

set.seed(12345678)
n = 100; beta0 = 2.5; beta1 = 0.8
x.lm = runif(n = n) # or it could be something else! rnorm(n = n, mean = 10, sd = 1)
err = rnorm(n = n, mean = 0, sd = 1)
# Linear combination 
mu = beta0 + beta1*x.lm 
y.lm = mu + err # same as y.lm = rnorm(n = n, mu, sd = 1)
# Make a dataframe of the data 
df.lm = data.frame(x = x.lm, y = y.lm)

round(coefficients(summary(lm.out)), 4)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   2.4726     0.1635 15.1185   0.0000
# x             0.8175     0.2779  2.9413   0.0041
# Adjusted R^2
summary(lm.out)$adj.r.squared
# [1] 0.0717415
# summary(lm.out)$fstatistic
# anova(lm.out)

GLM: Linear model

• You can simulate the response (fitted line), from a distribution used in a fitted model.

sim.lm = simulate(lm.out, nsim = 2000, seed = 12) # simulate based on predicted value and sigma of model
rx = range(c(rowMeans(sim.lm), fitted(lm.out)))
hist(rowMeans(sim.lm),xlim=rx,main="Hist simulation",xlab="");title(xlab="rowMeans(sim.lm)",line=2.2,cex.lab=1.2)
hist(fitted(lm.out), xlim = rx,main="Hist fitted",xlab="");title(xlab="fitted(lm.out)",line=2.2,cex.lab=1.2)

c(mean(rowMeans(sim.lm)), mean(y.lm),  mean(fitted(lm.out))) # compare
# [1] 2.892684 2.891622 2.891622
rbind(head(rowMeans(sim.lm)), head(fitted(lm.out))) # Average of all simulations (in row) gives the fitted values
#             1        2        3        4        5        6
# [1,] 3.182545 2.677072 3.248524 3.210164 2.470363 2.895410
# [2,] 3.187451 2.708652 3.246897 3.195649 2.479667 2.899208

GLM: Linear model

• Now that we have the core elements to generate a lm, we can put the simulation inside a function and vary some parameters
• If we are interested in how our confidence in the slope estimate changes based on sample size we could replicate the simulation a certain number of times
  • We could plot the density of the estimates that we'd get out of the simulations.
  • Below, slopes estimates were extracted from each simulation with various sample sizes.
  • As the sample size increases, our slope estimation gets closer to the "true" underlying value.

GLM: Linear model

• Same simulation as before, but now we plot the adjusted $R^2$ from the linear model.

• You can see from the plots that the effect size, the $R^2$, varies much more when the sample size is small. As the sample size increases, the spread around the effect size diminishes.

GLM: Linear model

• Simulate 2 independent variables ( $x_1$ and $x_2$ ).

# Coefficients table:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   5.4483     0.9409  5.7905        0
# x1            0.5419     0.0922  5.8803        0
# x2            0.4796     0.0157 30.5133        0
# R Squared:
# [1] 0.9083729

GLM: Linear model (with categorical)

• Simulate 1 numerical variable ( $x_1$ ) separated in 2 levels (Control, Treatment).

set.seed(12345678)
n = 100; meanx = 10
x.lm = rnorm(n = n, mean = meanx, sd = 1)
cat.lm = c("Control", "Treatment")
gr = gl(n = length(cat.lm), k = n/2, labels = cat.lm) # makes factors (default)
# Make a dummy matrix
dummy.mat = matrix(data = 0,
                   nrow = n,
                   ncol = length(cat.lm))
colnames(dummy.mat) <- paste0(c(cat.lm))
dummy.mat[,1] <- 1 # Intercept
dummy.mat[,2] <- ifelse(gr=="Treatment",1,0)
beta1 = 0.8; betaControl = .2; betaTreament = .7
# Combine the data to get the linear prediction 
lin.pred.x = dummy.mat %*% c(betaControl,betaTreament) + beta1 * x.lm
y.lm = lin.pred.x + rnorm(n = n, mean = 0, sd = 1)
df.lm = data.frame(x = x.lm, y = y.lm, gr = gr)
b.5 = scales::alpha("black",alpha = .5); r.5 = scales::alpha("red",alpha = .5)
lm.out = lm(y~x+gr, data = df.lm);   lmcoef <- lm.out$coefficients

GLM: Linear model (with categorical)

plot(y~x, data = df.lm, pch = 19, col = ifelse(gr=="Treatment",r.5,b.5))
abline(a= lmcoef["(Intercept)"], b= lmcoef["x"], col = "black") # control
abline(a= lmcoef["(Intercept)"] + lmcoef["grTreatment"], b=lmcoef["x"], col="red")
boxplot(y~gr, data = df.lm, pch = 19, col = ifelse(gr=="Treatment",r.5,b.5))

round(coefficients(summary(lm.out)), 4) # betaControl is (Intercept)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  -0.1184     1.0051 -0.1178   0.9065
# x             0.8449     0.0996  8.4855   0.0000
# grTreatment   0.7308     0.1861  3.9267   0.0002
summary(lm.out)$adj.r.squared           # betaTreament is grTreatment. 
# [1] 0.446509

GLM: Linear model (interaction; discrete)

• Linear models can also take interactions in consideration.
• $y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2+\epsilon$
• In this simple example $\beta_1$ and $\beta_2$ is the effect that the predictor $x_1$ and $x_2$ have, and $\beta_3$ is the effect that both $x_1$ and $x_2$ predictor have (or that the effect of one variable (e.g., $x_1$) depends on the value of the other (e.g., $x_2$)).

#             Estimate Std. Error  t value Pr(>|t|)
# (Intercept)   2.4321     0.2085  11.6661        0
# x1            1.5561     0.3094   5.0289        0
# x2            2.1088     0.1590  13.2626        0
# x1:x2        -3.0863     0.2276 -13.5584        0
# [1] 0.7190843

GLM: Linear model (interaction; discrete)

plot(x = df.lm[df.lm$x1 == 0, ]$x2, y = df.lm[df.lm$x1 == 0, ]$y, # Add x1 = 0
     col = rgb(red = 0, green = 0, blue = 1, alpha = 0.25), pch = 19,
     xlab = "x2", ylab = "y", ylim = range(df.lm$y)); abline(h=0, v=0, lt =3)
abline(a = coef(lm.out)[1], b = coef(lm.out)[3], col = "blue", pch = 19, lwd =2)
points(x = df.lm[df.lm$x1 == 1, ]$x2, y = df.lm[df.lm$x1 == 1, ]$y, # Add x1 = 1
       col = rgb(red = 1, green = 0, blue = 0, alpha = 0.25), pch = 19)
abline(a = coef(lm.out)[1] + coef(lm.out)[2], 
       b = coef(lm.out)[3] + coef(lm.out)[4], col = "red", lwd = 2)

GLM: Linear model (interaction; continuous)

• Here is the simulation of a linear model with an interaction term between two continuous variables.
• The plot is to show that the relationship between $y$ and $x_1$ is modulated by $x_2$
• If $\hat{y}=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2$, but you want to show the graph for $x_1$, we can rearrange the equation to this:
  • $\hat{y}=(\beta_0+\beta_2x_2)+(\beta_1+\beta_3x_2)x_1$
  • You can see from the equation above that in order to get an intercept and slope for $x_1$, you need to specify a $x_2$ value. Therefore, if you want to show the effect of $x_1$ on $y$ you need to set $x_2$ to a certain value. So for a given value of $x_2$, the intercept on the y-axis is the first term $(\beta_0+\beta_2x_2)$ while the slope of $x_1$ is $(\beta_1+\beta_3x_2)$.

GLM: Linear model (interaction; continuous)

hist(df.lm$x1);hist(df.lm$x2);hist(df.lm$y)

lm.out = lm(y~x1*x2, data = df.lm); lmcoef = coef(lm.out)
# Make a new range of x2 values on which we will test the effect of x1 
x2r = range(x.lm2); x2.sim = seq(x2r[1], x2r[2], by = .5); x2s.lg=length(x2.sim)
# This is the effect of x1 at different values of x2 (moderates the effect of x1)
eff.x1 <- lmcoef["x1"] + lmcoef["x1:x2"] * x2.sim # gets slopes  
eff.x1.int <- lmcoef["(Intercept)"] + lmcoef["x2"] * x2.sim # gets  int.  
eff.dat <- data.frame(x2.sim, eff.x1, eff.x1.int) # Put that in dataframe 
virPal <- viridis::viridis(length(x2.sim), alpha = .8) # Get colours 
eff.dat$col <- virPal[as.numeric(cut(eff.dat$x2.sim, breaks = x2s.lg))]
df.lm$col <- virPal[as.numeric(cut(df.lm$x2, breaks = x2s.lg))]
df.lm$line <- c("black","red")[as.numeric(cut(df.lm$x2, breaks = x2s.lg)) %% 2+1]

GLM: Linear model (interaction; continuous)

par(mfrow=c(1,1), mar =c(4,4,1,1))
plot(x = df.lm$x1, y = df.lm$y, bg = df.lm$col, 
     pch = 21, xlab = "x1", ylab = "y")
apply(eff.dat, 1, function(x) abline(a = x[3], b = x[2], col = x[4], lwd  = 2))
# NULL
abline(h = 0, v = 0,lty = 3)
legend("topleft", title = "x2",legend = round(eff.dat$x2.sim,1), lty = 1, lwd = 3,
       col = eff.dat$col, bg = scales::alpha("white",.5))

Challenge 6

• You have 20 seeds from plants that are genetically similar on which you want to see the effect of a treatment, say a species of fungi added to the soil compared to a control with only the soil.
• After a certain time, say 6 months, you want to weight the dried plant material and test if there is an effect of your treatment
• You've read in the literature that similar experiments have had a mean weight of 5.03g (standard deviation = 0.28) for the control and 4.66g (standard deviation = 0.49) for the treatment.
• Your scale adds and error of 0.3g of standard deviation to the measurement.
• Can you tell the proportion of p-values that will be significant with these details?

Challenge 6 - Solution

set.seed(1234)
n = 10;  sd.e = 0.3
control = rnorm(n, mean = 5.03, sd = 0.28)
treatmt = rnorm(n, mean = 4.66, sd = 0.49)
gr <- gl(n = 2, k = n, length = n*2, labels = c("ctl","trt"))
weight = c(control,treatmt) + rnorm(2*n, mean = 0, sd = sd.e)
plant.weight = data.frame(weight=weight, gr = gr)
lm.out.int <- lm(weight ~ gr, data = plant.weight)
lm.out.noi <- lm(weight ~ gr-1, data = plant.weight) # Comparing the data to 0 (int. = 0, which becomes the reference)
anova(lm.out.int) # this time, it was significant 
# Analysis of Variance Table
# 
# Response: weight
#           Df Sum Sq Mean Sq F value  Pr(>F)  
# gr         1 0.9422 0.94219  4.8569 0.04079 *
# Residuals 18 3.4918 0.19399                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summary(lm.out.noi) # Just to see the estimate (NOT THE P-VALUES)
s.lm.out = summary(lm.out.int)
# You can get the p-value with this function 
pt(q = abs(s.lm.out$coefficients["grtrt","t value"]),df = s.lm.out$df[2], lower.tail = F) *2
# [1] 0.04079444

Challenge 6 - Solution

plot(weight~gr, data=plant.weight, las=1)

plot(lm.out.int, las = 1)

Challenge 6 - Solution

# Make the function for the simulation
exp.plant <- function(n = 100, # number of seeds in EACH group (2*n = total)
                      sd.e = 0.3, plot=F, ret = T) {
  control = rnorm(n, mean = 5.03, sd = 0.28)
  treatmt = rnorm(n, mean = 4.66, sd = 0.49)
  gr <- gl(n = 2, k = n, length = 2*n, labels = c("ctl","trt"))
  weight = c(control,treatmt) + rnorm(2*n, mean = 0, sd = sd.e)
  plant.weight = data.frame(weight=weight, gr = gr)
  if (plot) {
    plot(weight ~ gr, data = plant.weight, las = 1)
  }
  lm.out.int <- lm(weight ~ gr, data = plant.weight)
  s.lm.out = summary(lm.out.int)
  if (ret) {
   return(s.lm.out) 
  }
}

Challenge 6 - Solution

set.seed(2345); nb.rep = 2000 # number of replications in the simulation 
l.rp = replicate(n = nb.rep, simplify = FALSE,
                 expr = exp.plant( n = 10)$coefficients["grtrt","t value"])
p.val.lm = pt(q = abs(unlist(l.rp)),df = s.lm.out$df[2], lower.tail = F)*2 # Get p-values 
exp.plant(n = 10, plot = T, ret = F) # plot a simulation example 
hist(unlist(l.rp), main = "t-values", probability = T, xlab ="")
lines(density(unlist(l.rp)), col="blue", lwd=2); abline(v = qt(0.025, df = s.lm.out$df[2]))
hist(p.val.lm, main = "p-values", probability = T, xlab ="", xlim = c(0,1))
lines(density(p.val.lm), col="red", lwd=2); abline(v = 0.05)

• The proportion of p-values that are significant is 31.45%

GLM: Logistic

• Recall that logistic is $Y \sim Binomial(p)$ with $\text{log} \Bigl(\frac{p}{1-p}\Bigr) = \beta_0 +\beta_1x$ or $p=\frac{e^{\beta_0+\beta_1x}}{1+e^{\beta_0+\beta_1x}}$ where $\text{log} \Bigl(\frac{p}{1-p}\Bigr)$ is the log odds or log likelihood. The Y values are determined by a Bernoulli distribution (binomial of size = 1)
• We use maximum likelihood to fit a line to the data. It doesn't have 'residuals' and don't have an $R^2$
• Remember that in order to make the domain $]-\infty,\infty[$ for the y values (instead of $[0,1]$), we need to transform the probability to the $log(odds)$.

GLM: Logistic

• This simulation concerns a continuous variable.

set.seed(987654)
n = 100
x1 = rnorm(n = n, mean = 6, sd = 1)
x2 = rnorm(n = n, mean = 0, sd = 1)
# Rescale the data
x1z = scale(x1)
x2z = scale(x2)
z = 0 + 2*x1z + 3*x2 # This gives the LOG odds. Recall that $p = odds/(1+odds)$
pr = 1/(1+exp(-z)) # Transform to get the LOG(odds) # inverse-logit function; 
y = rbinom(n = n, size = 1, prob = pr) # Bernoulli response variable 
# Combine the data in a dataframe 
df = data.frame(y = y, x1 = x1, x2 = x2)
# Compute the model
glm.logist = glm( y~x1+x2, data=df, family="binomial")

GLM: Logistic

plot(y~x1, data = df, col = scales::alpha("black",.5), pch = 19)
newdata <- data.frame(x1=seq(min(x1), max(x1),len=n), 
                      x2 = seq(min(x2), max(x2),len=n))
newdata$y = predict(object = glm.logist, newdata = newdata, type = "response") 
lines(x = newdata$x1, y = newdata$y, col = "red",lwd = 2)

• The coefficients are given in log(odds)

# (Intercept)          x1          x2 
#  -11.503498    1.987659    3.061898

• The z-values are the number of sd the estimate is from a z-distribution.
• For every one unit in X, the log(odds) increases by about 2.

GLM: Logistic

• Here is another example of the logistic regression with only 1 independent variable (z = 0 + 2*rnorm(n = 20, mean = 6, sd = 1))
• You can see that the coefficients are given in log odds.

#              Estimate Std. Error   z value     Pr(>|z|)
# (Intercept) -4.752291  0.7031429 -6.758641 1.392917e-11
# x1           0.835486  0.1179451  7.083686 1.403694e-12

# Probability at intercept is 0.008558029
# Probability for an increase in 1 sd of X is 0.6975137

GLM: Logistic (quadratic)

#               Estimate Std. Error   z value     Pr(>|z|)
# (Intercept)  1.1567584  0.2365723  4.889662 1.010094e-06
# x1           0.2510816  0.1780770  1.409961 1.585512e-01
# I(x1^2)     -0.3458782  0.1378428 -2.509222 1.209974e-02
# [1] 1
# Probability at intercept is 0.3339322

GLM: Logistic (categories)

• We can also simulate a model with categorical variables

GLM: Poisson

• Recall that Poisson is $Y \sim Poisson(\mu)$ with $\text {ln} \mu=\beta_0+\beta_1x$ or $\mu=e^{\beta_0+\beta_1x}$ (no separate error term as $\lambda$ determines both the mean and variance)

set.seed(42)
n = 500
x = rnorm(n = n, mean = 0, sd = 1)
# Rescale the data
xz = scale(x)
log.mu = 1 + .8*xz
y = rpois(n = n, lambda = exp(log.mu)) 
# Combine the data in a dataframe 
df = data.frame(y = y, x = x)

GLM: Poisson

#now feed it to glm:
glm.poisson = glm( y~x, data=df, family="poisson")
plot(y~x, data = df, col = scales::alpha("black", 0.5), pch = 19, cex = 0.5,
     main = "Poisson regression")
newdata <- data.frame(x = seq(min(x), max(x), len = n))
newdata$y = predict(object = glm.poisson, newdata = newdata, type = "response") 
lines(x = newdata$x,
      y = newdata$y, col = "red",lwd = 2)
glm.gau = glm( y~x, data=df, family="gaussian")
hist(residuals(glm.gau), main = "Residuals Gaussian")
hist(residuals(glm.poisson), main = "Residuals Poisson")

Linear mixed model (LMMs) refresher

• A simple linear model is actually the 'simplest' mixed model (although the convention is that we don't called it mixed)
• What is the random effect in a linear model?

$$Y = \beta_{0} + \beta_{1} x_{1} + \cdots + \beta_{p} x_{p} + \epsilon$$

Linear mixed model (LMMs) refresher

• A simple linear model is actually the 'simplest' mixed model (although the convention is that we don't called it mixed)
• What is the random effect in a linear model?

$$Y = \beta_{0} + \beta_{1} x_{1} + \cdots + \beta_{p} x_{p} + \epsilon$$

• The residual of the model ( $\epsilon$ ) is actually the 'random effect' since it is drawn from a distribution which can change for each data point $\epsilon \sim N(\mu=0,sd = \sigma)$.

Linear mixed model (LMMs) refresher

• to simplify the linear model equation we are going to write it this way: $$\mathbf{Y} = \mathbf{X}\beta + \epsilon$$

• This only means that $\mathbf{Y}$ is a vector (or matrix) made up of the linear combination of fixed effects $\mathbf{X}\beta$ and a random part $\epsilon$.

• In this type of notation, LMMs is (random intercept): $$\mathbf{Y} = \mathbf{X}\beta + \mathbf{Z\upsilon} + \epsilon$$

• With all the previous elements being the same but the random part which is composed of $\mathbf{Z \upsilon} + \epsilon$, for all categories $\mathbf{\upsilon}$. $\mathbf{X}$ is the design matrix and $\mathbf{Z}$ is the block matrix.

• Just keep in mind that, as before, $\epsilon \sim N(\mu=0,sd = \sigma)$ and $\mathbf{\upsilon} \sim N(\mu=0,sd = \mathbf{D})$, where $\mathbf{D}$ is a variance-covariance matrix. Also, $\mathbf{\upsilon}$ and $\epsilon$ are independent.

Linear mixed model (LMMs) refresher

• This is just another way to write the model (there is a part that specify the random slopes)

$$y_{ij} = \beta_{0} + \beta_{1} x_{1ij} + \upsilon_{1j} x_{1ij} + \upsilon_{0ij} + \epsilon_{0ij}$$

• $\epsilon_{0ij} \sim N(\mu=0,sd = \sigma_{e0})$ and $\mathbf{\upsilon_{0j}} \text{ and } \mathbf{\upsilon_{1j}} \sim N(\mu=0,sd = \Omega_{\upsilon})$ and $\Omega_{\upsilon} = \left[\begin{align*} \sigma_{\upsilon0}^2 & \sigma_{\upsilon10} \\ \sigma_{\upsilon01} & \sigma_{\upsilon1}^2 \end{align*}\right]$

Linear mixed model (LMMs) refresher

• Just so that everything is extra clear $$\mathbf{Y} = \mathbf{X}\beta + \mathbf{Z \upsilon} + \epsilon$$

One implementation of that could be ( $\beta_1$ and $\beta_2$ would represent different treatments, $v_1$ and $v_2$ are different categories):

$$\left[\begin{array}{l} y_{11} \\ y_{21} \\ y_{12} \\ y_{22} \end{array} \right] = \left[ \begin{array}{ll} 1 & 0 \\ 0 & 1 \\ 1 & 0 \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} \beta_{1} \\ \beta_{2} \end{array}\right]+\left[\begin{array}{ll} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} v_{1} \\ v_{2} \end{array}\right]+\left[\begin{array}{l} \epsilon_{11} \\ \epsilon_{21} \\ \epsilon_{12} \\ \epsilon_{22} \end{array}\right]$$

LMMs simulation

• (G)LMMs are a neat extention of the LM models that can take into account different type of response variables and variance parterns in the response data (see QCBS workshop on the topic)