1.1.1 Deterministic model

• Generate 100 samples of \(x,y\) by
  • generate genotypes \(x\) for each sample from a uniform distribution on integers \(0,1,2\) (tip: use runif)
  • set \(\beta_0=0.3\) and \(\beta_1=0.2\)
  • generate a phenotype \(y\) for each sample of using the model \(y = \beta_0 + \beta_1* x\)
  • plot \(y\) against \(x\).

#parameters
b0 = 0.3
b1 = 0.2
N=100
sim1 = data.frame(x=round(runif(N,min=0,max=2)))

sim1$y = b0 + b1 * sim1$x

plot(x=sim1$x, y=sim1$y)

1.1.2 Statistical model

Deterministic vs statistical models

• \(Y = \beta_0 + \beta_1* X\) is a deterministic model
  • does not model any variation
  • Common, e.g., in classical physics (velocity \(v = \frac{\Delta x}{\Delta t}\))
• \(Y = \beta_0 + \beta_1* X +\epsilon,\) where \(\epsilon~\sim N(mean =0,sd=\sigma)\) is a statistical (equiv. stochastic, random) linear model (\(\sim\) means ‘is drawn from’)
  • attempts to model variation around a population mean (the residuals) determined by the model
    • the residual of sample \(i\) is \(y_i-(\beta_0 + \beta_1* X)\).
  • used in statistic analysis
  • \(Y\) is then called a stochastic or random variable

Henceforth, we will only consider statistical models.

Task

• Generate data from \(Y = \beta_0+ \beta_1* X +\epsilon,\) where \(\epsilon ~\sim N(mean=0,sd = 0.05)\) (tip: use rnorm).

#parameters
b0 = 0.3
b1 = 0.2
N=100
sim1 = data.frame(x=round(runif(N,min=0,max=2)))

sim1$y = b0 + b1 * sim1$x + rnorm(N, mean=0, sd=0.05) 

plot(x=sim1$x, y=sim1$y)

1.2.1 Simulation solution

• simulate 1000 \(y\) using the model above with \(x=2\) and store in a vector \(Y\).
• Make a histogram of \(Y\) and estimate the probability that \(Y<=0.65\) (tip: make sure that \(0.65\) is among to a histogram breakpoints)

Think About

• What shape does the plotted histogram have?
• Where approximately is the mean?
• Does this make sense in light of the generative model we used?
• Can we plot \(Pr[Y<=y|X,\theta]\) more directly? (tip: google, e.g., ‘probability distribution’)

#parameters
b0 = 0.3
b1 = 0.2
N = 1000
x = 2

y = b0 + b1 * x + rnorm(N, mean=0, sd=0.05) 

h=hist(y, frequency=TRUE, right=FALSE, breaks=seq(0,1.0, 0.05), plot=FALSE)
h$counts = h$counts/N
plot(h, xlim=c(0,1.0), ylim =c(0,1), labels=TRUE)

# compute requested probability from hist
paste("Pr[Y<=0.65|X=2] = ", sum(h$counts[1:13]))

#parameters
b0 = 0.3
b1 = 0.2
N = 1000
x = 2

y = b0 + b1 * x + rnorm(N, mean=0, sd=0.05) 

h=hist(y, frequency=TRUE, right=FALSE, breaks=seq(0,1.0, 0.05), plot=FALSE)
h$counts = h$counts/N

# replace the histogram counts with the cumulative counts
h$counts = cumsum(h$counts)
#plot the CDF
plot(h, labels=TRUE)

# compute requested probability from hist
paste("Pr[Y<=0.65|X=2] = ", sum(h$counts[13]))

## [1] "Pr[Y<=0.65|X=2] =  0.151"

1.2.2 Analytical solution

• Conveniently, R has both the PDF and CDF implemented, making a
• Calculate \(Pr[Y<=0.65|X=2]\) using R’s cumulative distribution function (CDF) for the Normal distribution (tip: look up help for Normal)
  • The CDF of the value \(y\) is the probability that the stochastic variable \(Y<=y\).
• Challenge: try both approaches with other values of \(y\) or try to compute \(Pr[Y>y|X=2]\).

Think about

• Do the result fit that from the Simulation?
  • If not, can make the fit better?
• What can this result be used for?
• Challenge: Why do we look at intervals of \(Y\) rather than specific values?

b0 = 0.3
b1 = 0.2
x = 2
y = 0.65
mu = b0 + b1 * x

paste("Pr[Y<=.6|X=2] = ",pnorm(y, mean = mu, sd=0.05))

1.3.1 Student’s t-test

We will use a single sample t-test, which tests how probable it is that some sampled data is drawn from a given Normal distribution; this distribution is called the NULL model of the test. The t-test is not designed to handle different genotypes, so for simplicity, we will, again, limit ourselves to the case \(X=2\).

The t-test uses normalized residuals as a test statistics:

\[\begin{eqnarray*} t &= \sum_{i}\frac{y_i-\mu}{s/\sqrt{N}} &= \frac{\bar{y}-\mu}{s/\sqrt{N}} \end{eqnarray*}\] where \(s/\sqrt{N}\) is an estimate of the standard deviation, \(\sigma\), from the observed data. As you can see this boils down to comparing the (standardized) estimated mean from the data, \(\bar{y}\) and the model mean, \(\mu\).

• Challenge: show the equality above`

These normalized residuals are approximately distributed \(t~\sim N(0,1)\) (this is the General Normal distribution, see also ‘Possible answers/Extra Reading’ under 1.2.1)

Conveniently, The t-test is implemented in the R function t.test

• Simulate a new data set from a linear model \(Y' = \gamma_0 + \gamma_1 X'\), for \(X'=2\) for all samples.
  • For simplicity, we will limit ourselves to the case \(X'=2\)
• Use our previous model \(Y=\beta_0+\beta_1 X\) as the NULL model.
  • Perform a one sample t-test (with \(\mu\) from the NULL model) on new \(Y'\) data.
• Try some different values of \((\gamma_0, \gamma_1)\) closer or further away from \((\beta_0, \beta_1)\)

Think about

• Does the t-test reject the NULL model or not?

# compute the mean, mu, from our model 
b0 = 0.3
b1 = 0.2

mu = b0 + b1 * x

# Let's simulate the data from another model
c0 = 0.4
c1 = 0.2
N=100
yp = c0 + c1 * rep(x,N) + rnorm(N, mean=0, sd=0.05) 
# t-test
t.test(yp, mu=mu)

# data 1
b0 = 0.3
b1 = 0.2
N=100
y = b0 + b1 * rep(x,N) + rnorm(N, mean=0, sd=0.05) 

#data 2
c0 = 0.4
c1 = 0.2
N=100
yp = c0 + c1 * rep(x,N) + rnorm(N, mean=0, sd=0.05) 

t.test(y,yp)

# reuse previously simulated data in sim1
b0 = 0.3
b1 = 0.2
N=100
sim1 = data.frame(x=round(runif(N,min=0,max=2)))
sim1$y = b0 + b1 * sim1$x

summary(lm(y~x, data=sim1))