Chapter 4 - Geocentric Models

This chapter introduced the simple linear regression model, a framework for estimating the association between a predictor variable and an outcome variable. The Gaussian distribution comprises the likelihood in such models, because it counts up the relative numbers of ways different combinations of means and standard deviations can produce an observation. To fit these models to data, the chapter introduced quadratic approximation of the posterior distribution and the tool quap. It also introduced new procedures for visualizing prior and posterior distributions.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.

Questions

4-1. For the model definition below, simulate observed y values from the prior (not the posterior). Make sure to plot the simulation. \[\begin{align} \ y_i ∼ Normal(μ, σ) \\ \ μ ∼ Normal(0, 10) \\ \ σ ∼ Exponential(1) \\ \end{align}\]

sample_mu <- rnorm(1e4, 0, 10)
sample_sigma <- runif(1e4, 0, 10)
prior_y <- rnorm(1e4, sample_mu, sample_sigma)
dens(prior_y)

4-2. Translate the model just above into a quap formula.

formula <- alist(
  y ~ dnorm(mu, sigma),
  mu ~ dnorm(0, 10),
  sigma ~ dunif(0, 10)
)
4-3. Translate the quap model formula below into a mathematical model definition:

y ~ dnorm( mu , sigma ),
mu <- a + b*x,
a ~ dnorm( 0 , 10 ),
b ~ dunif( 0 , 1 ),
sigma ~ dexp( 1 )

flist <- alist(
  y ~ dnorm(mu, sigma),
  mu <- a + b*x,
  a ~ dnorm(0, 50),
  b ~ dunif(0, 10),
  sigma ~ dunif(0, 50)
)

4-4. A sample of students is measured for height each year for 3 years. After the third year, you want to fit a linear regression predicting height using year as a predictor. Write down the mathematical model definition for this regression, using any variable names and priors you choose. Be prepared to defend your choice of priors. Simulate from the priors that you chose to see what the model expects before it sees the data. Do this by sampling from the priors. Then consider 50 students, each simulated with a different prior draw. For each student simulate 3 years. Plot the 50 linear relationships with height(cm) on the y-axis and year on the x-axis. What can we do to make these priors more likely?

# A linear model is chosen. A normal distribution with an SD of 25 cm for the alpha prior; 150 cm is in the middle of the expected distribution if both high school and college students are included, and 25 cm is enough variability that two SDs around the mean (i.e., 100 cm to 200 cm) should include most students at the high and low end of the age distribution.
# A normal distribution with SD of 2 cm/year for the beta prior; 4 cm/year is in the middle of the expected distribution if both high school and college students are included, and 2 cm/year is enough variability that two SDs around the mean (i.e., 0 cm/year to 8 cm/year) should include most students at the high and low end of the age distribution. 
#And for the sigma before, A uniform distribution ranging from 0 to 50 cm; this range encompasses both a close distribution of kids around the same age/height as well as a large range of students at both school and college ages/heights.
# hi∼Normal(μ,σ)    
# μi = α+βxi
# α ∼Normal(150,25)
# β∼Normal(4,2)
# σ ∼Uniform(0,50)
dens(rnorm(1e4, 150, 25))

dens(rnorm(1e4, 4, 2)) 

dens(rnorm(1e4, 150, 30))

4-5. (a) Now suppose I remind you that every student got taller each year. Does this information lead you to change your choice of priors? How? Again, simulate from the priors and plot. (b) Now suppose I tell you that the variance among heights for students of the same age is never more than 64cm. How does this lead you to revise your priors?

# This information provides more confidence that we're dealing with school-aged children (e.g., around 7 years old). As a consequence, we may narrow the range of  past distributions to minimize the likelihood of older ages' heights and growth rates. Our expectations  are also considerably lower now, as an equal mix of junior and senior students may not be expected.

# hi∼Normal(μ,σ)    
# μi = α+βxi
# α ∼Normal(120,10)
# β∼Normal(7,1)
# σ ∼Uniform(0,20)



N <- 1000
a <- rnorm( 1e4, 120, 10)
b <- rnorm( N, 7, 1 )
sigma <- runif( N, 0, 20 )

height_1 <- rnorm(N, a + b * 0, sigma)
dens(height_1)

4-6. In the chapter, we used 15 knots with the cherry blossom spline. Increase the number of knots and observe what happens to the resulting spline. Then adjust also the width of the prior on the weights—change the standard deviation of the prior and watch what happens. What do you think the combination of knot number and the prior on the weights controls?

data(cherry_blossoms) 
d <- cherry_blossoms 
d2 <- d[ complete.cases(d$temp) , ]
# Increase the number of knots
numknots <- 30
knotlist <- quantile( d2$year , probs=seq(0,1,length.out=numknots))
library(splines) 
B <- bs(d2$year, knots=knotlist[-c(1,numknots)] , degree=3 , intercept=TRUE )
plot( NULL , xlim=range(d2$year) , ylim=c(0,1) , xlab="year" , ylab="value" )
for ( i in 1:ncol(B) ) lines( d2$year , B[,i])

# The curve of 30 knots is wigglier than the curve with 15 knots.