I’m trying to draw a regression line like this picture . However, the result of my code looks like this.
What am I missing or where am I wrong? Thank you so much!
data("cherry_blossoms")
d <- cherry_blossoms
d2 <- d[complete.cases(d), ]
num_knots <- 15
knot_list <- quantile( d2$temp, probs = seq(0, 1, length.out = num_knots ) )
B <- bs(d2$temp, knots=knot_list[-c(1, num_knots)],degree=3, intercept=TRUE)
spline_mdl <- quap(
alist(
D ~ dnorm( mu, sigma ),
mu <- a + B %*% w,
a ~ dnorm( 100, 10 ),
w ~ dnorm( 0, 10 ),
sigma ~ dexp(1)
),
data = list( D = d2$doy , B = B) ,
start = list( w=rep( 0 , ncol(B)))
)
mu <- link( spline_mdl , data = d2)
mu_mean <- apply ( mu , 2 , mean )
mu_PI <- apply ( mu , 2 , PI , prob = 0.89)
doy_sim <- sim ( spline_mdl , data = d2)
doy_PI <- apply( doy_sim , 2 , PI , prob = 0.89)
plot(d2$temp , d2$doy , col = col.alpha(rangi2, 0.5), xlab = "Temperature",
ylab = "Day of Year" )
lines(d2$temp , mu_mean)
shade( mu_PI , d2$temp)
shade( doy_PI , d2$temp)
mtext("Splines with 15 knots")
The question I’m practicing is the following:
Model the association between blossom date (doy) and March temperature (temp). Note that there are many missing values in both variables. You may consider a linear model, a polynomial, or a spline on temperature. How well does temperature trend predict the blossom trend?
- Operating System: macOS