We fit a linear regression model to the !Kung data within the Rethinking package using 6 methods. 1) R’s builtin lm(), 2) likelihood function, 3) Rethinking quap, 4) Rethinking ulam, 5) Nimble, and 6) Rstan/Stan.

This exercise uses the !Kung height vs. weight data from the Rethinking package for adults age 18 or older.

Load and examine the data.

library(rethinking)
## Loading required package: rstan
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Loading required package: parallel
## rethinking (Version 2.13)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:stats':
## 
##     rstudent
data(Howell1)
d <- Howell1
d2 <- d[ d$age >= 18 , ]
head(d2)
##    height   weight age male
## 1 151.765 47.82561  63    1
## 2 139.700 36.48581  63    0
## 3 136.525 31.86484  65    0
## 4 156.845 53.04191  41    1
## 5 145.415 41.27687  51    0
## 6 163.830 62.99259  35    1
plot(height~weight, data=d2)

# Scaling the data to mean = 0, sd = 1.
d3<-d2
d3[c(1,2)] <- lapply(d2[c(1, 2)], function(x) c(scale(x)))

# Equivalent code but more transparent
d4<-d2
d4[1]<-scale(d4[1])
d4[2]<-scale(d4[2])

# identical(d4[,1],d3[,1])
# all.equal(d4[,1],d3[,1])
# sum(d3[,1]-d3[c(1)])
#sum(d3[,1]-d4[c(1)])
# all.equal(d3[,1],d4[,1])
  1. Fit a linear regression to these data using the ‘lm’ command within R. The model should predict the observed weight from height and the model should include an intercept and an effect of height. Plot the estimated linear fit superimposed on the data and report the parameter estimates for the intercept, effect of height, and sigma.
out<-lm(height~weight,data=d2)
summary(out,correlation=TRUE) # # correlation about -0.99
## 
## Call:
## lm(formula = height ~ weight, data = d2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.7464  -2.8835   0.0222   3.1424  14.7744 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 113.87939    1.91107   59.59   <2e-16 ***
## weight        0.90503    0.04205   21.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.086 on 350 degrees of freedom
## Multiple R-squared:  0.5696, Adjusted R-squared:  0.5684 
## F-statistic: 463.3 on 1 and 350 DF,  p-value: < 2.2e-16
## 
## Correlation of Coefficients:
##        (Intercept)
## weight -0.99
# Repeating but with scaled columns
out2<-lm(height~weight,data=d3)
summary(out2,correlation=TRUE) # correlation now 0
## 
## Call:
## lm(formula = height ~ weight, data = d3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.55045 -0.37244  0.00287  0.40587  1.90826 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.794e-17  3.502e-02    0.00        1    
## weight      7.547e-01  3.507e-02   21.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.657 on 350 degrees of freedom
## Multiple R-squared:  0.5696, Adjusted R-squared:  0.5684 
## F-statistic: 463.3 on 1 and 350 DF,  p-value: < 2.2e-16
## 
## Correlation of Coefficients:
##        (Intercept)
## weight 0.00
plot(height~weight, data=d3)
abline(a=out2$coefficients[1],b=out2$coefficients[2])

library(stats)
sigma(out2)
## [1] 0.6569515
  1. We fit the linear model using our own likelihood function rather than R’s built in ‘lm’. Plot the estimated linear fit superimposed on the data and report the parameter estimates for the intercept, effect of height, and sigma. Are they the same as from ‘lm’?

Here is our likelihood function. We actually work with the negative log likelihood. Why?

nllNormReg<-function(parVec,data=d3){
  wt<-data$weight
  ht<-data$height
  intercept<-parVec[1]
  beta<-parVec[2]
  mySd<-parVec[3]
  htPred<-intercept+beta*wt
  nllik<- -sum(dnorm(x=ht,mean=htPred,sd=mySd,log=TRUE))
  #browser()
  return(nllik)
}

Minimizing the negative log likelihood function above. This is equivalent to finding the maximum likelihood estimate.

parVec<-c(1.0,2.0,2.0) # Need some starting parameters for gradient climbing/descending methods
outLM<-optim(par=parVec,fn=nllNormReg,method="L-BFGS-B",
               lower=c(-Inf,-Inf,0.0001),upper=c(Inf,Inf,Inf))
outLM$par # 
## [1] -2.940787e-07  7.547480e-01  6.550838e-01

From R’s lm above: intercept = 5.794e-17 (~ 0) weight = 7.547e-01 sigma = 0.6569515

  1. We fit the linear model using the quap in the Rethinking package. Report the parameter estimates for the intercept, effect of height, and sigma. Are they the same as from ‘lm’? The quap function uses a normal approximation of posterior.
library(rethinking)

modelQUAP <- quap(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- intercept + beta*weight ,
        intercept ~ dnorm( 10 , 20 ) ,
        beta ~ dnorm( 0 , 5 ) ,
        sigma ~ dunif( 0 , 50 )
    ) , data=d3 )

precis( modelQUAP )
##                   mean         sd        5.5%      94.5%
## intercept 3.056654e-05 0.03491596 -0.05577189 0.05583302
## beta      7.547110e-01 0.03496486  0.69883036 0.81059157
## sigma     6.550825e-01 0.02468910  0.61562459 0.69454048

From R’s lm above: intercept = 5.794e-17 (~ 0) weight = 7.547e-01 sigma = 0.6569515

  1. We fit the linear model using the ulam function in the Rethinking package. Report the parameter estimates for the intercept, effect of height, and sigma. Are they the same as from ‘lm’? The ulam function uses MCMC sampling to estimate the posterior.
library(rethinking)

modelMCMC <- ulam(
      alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- intercept + beta*weight ,
        intercept ~ dnorm( 178 , 20 ) ,
        beta ~ dnorm( 0 , 1 ) ,
        sigma ~ dunif( 0 , 50 )
      )
   , data=list(height=d3$height, weight=d3$weight) , chains=3 )
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '1dd7f4e218caedbf89a4575a52ebc4c8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.045894 seconds (Warm-up)
## Chain 1:                0.041811 seconds (Sampling)
## Chain 1:                0.087705 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '1dd7f4e218caedbf89a4575a52ebc4c8' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.047635 seconds (Warm-up)
## Chain 2:                0.041745 seconds (Sampling)
## Chain 2:                0.08938 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '1dd7f4e218caedbf89a4575a52ebc4c8' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.042523 seconds (Warm-up)
## Chain 3:                0.041916 seconds (Sampling)
## Chain 3:                0.084439 seconds (Total)
## Chain 3:
precis(modelMCMC)
##                   mean         sd        5.5%      94.5%    n_eff     Rhat4
## intercept -0.001911767 0.03441552 -0.05915048 0.05383798 1238.072 0.9985964
## beta       0.753505680 0.03459942  0.69674576 0.80712866 1331.061 0.9996512
## sigma      0.659191367 0.02488534  0.62016521 0.69912745 1638.452 0.9996442

From R’s lm above: intercept = 5.794e-17 (~ 0) weight = 7.547e-01 sigma = 0.6569515

  1. We fit the linear model using the Nimble package in R. Report the parameter estimates for the intercept, effect of height, and sigma. Are they the same as from ‘lm’?
library(nimble)
## nimble version 0.11.1 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.
## 
## Attaching package: 'nimble'
## The following objects are masked from 'package:rethinking':
## 
##     cloglog, logit
## The following object is masked from 'package:stats':
## 
##     simulate
myCode <- nimbleCode({
  intercept ~ dnorm(0, sd = 10)
  beta ~ dnorm(0, sd = 10)
  sigma ~ dunif(0, 50)  # prior for variance components based on Gelman (2006)
  for(i in 1:n) {
    mu[i]<-intercept + beta*weight[i]
    height[i] ~ dnorm(mu[i], sd = sigma) 
  }
})

myConstants <- list(n = dim(d3)[1])
myData <- list(weight = d3$weight, height = d3$height)
myInits <- list(intercept = mean(d3$weight), beta = 0, sigma = sqrt(var(d3$weight)))
mcmc.out <- nimbleMCMC(code = myCode, data = myData, inits = myInits,
                       monitors = c("intercept", "beta","sigma"), 
                       thin = 10,
                       niter = 10000, nburnin=1000,
                       nchains = 1, setSeed = FALSE,
                       constants = myConstants,samplesAsCodaMCMC=TRUE)
## defining model...
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## checking model calculations...
## model building finished.
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
summary(mcmc.out)
## 
## Iterations = 1:900
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 900 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                 Mean      SD  Naive SE Time-series SE
## beta       0.7541926 0.03440 0.0011468      0.0011468
## intercept -0.0003289 0.03496 0.0011653      0.0011653
## sigma      0.6601961 0.02629 0.0008762      0.0008999
## 
## 2. Quantiles for each variable:
## 
##               2.5%      25%      50%     75%   97.5%
## beta       0.68936  0.73003 0.754140 0.77774 0.82236
## intercept -0.07228 -0.02278 0.001988 0.02256 0.06499
## sigma      0.61195  0.64184 0.658802 0.67791 0.71561

From R’s lm above: intercept = 5.794e-17 (~ 0) weight = 7.547e-01 sigma = 0.6569515

library(coda)
## 
## Attaching package: 'coda'
## The following object is masked from 'package:rstan':
## 
##     traceplot
traceplot(mcmc.out)

densplot(mcmc.out)

  1. We fit the linear model using the rstan package in R, which is an interface to Stan (outside of R). Report the parameter estimates for the intercept, effect of height, and sigma. Are they the same as from ‘lm’?
library("rstan")
# options(mc.cores = parallel::detectCores())
# rstan_options(auto_write = TRUE)

kung_dat <- list(N=dim(d3)[1],weight = d3$weight,height = d3$height)

fit <- stan(file = 'lm.stan', data = kung_dat,
            iter=10000, warmup=2500,thin=10,chains=4)
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL 'lm' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.53 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 2501 / 10000 [ 25%]  (Sampling)
## Chain 1: Iteration: 3500 / 10000 [ 35%]  (Sampling)
## Chain 1: Iteration: 4500 / 10000 [ 45%]  (Sampling)
## Chain 1: Iteration: 5500 / 10000 [ 55%]  (Sampling)
## Chain 1: Iteration: 6500 / 10000 [ 65%]  (Sampling)
## Chain 1: Iteration: 7500 / 10000 [ 75%]  (Sampling)
## Chain 1: Iteration: 8500 / 10000 [ 85%]  (Sampling)
## Chain 1: Iteration: 9500 / 10000 [ 95%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.161867 seconds (Warm-up)
## Chain 1:                0.552795 seconds (Sampling)
## Chain 1:                0.714662 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'lm' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 2501 / 10000 [ 25%]  (Sampling)
## Chain 2: Iteration: 3500 / 10000 [ 35%]  (Sampling)
## Chain 2: Iteration: 4500 / 10000 [ 45%]  (Sampling)
## Chain 2: Iteration: 5500 / 10000 [ 55%]  (Sampling)
## Chain 2: Iteration: 6500 / 10000 [ 65%]  (Sampling)
## Chain 2: Iteration: 7500 / 10000 [ 75%]  (Sampling)
## Chain 2: Iteration: 8500 / 10000 [ 85%]  (Sampling)
## Chain 2: Iteration: 9500 / 10000 [ 95%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.159085 seconds (Warm-up)
## Chain 2:                0.621557 seconds (Sampling)
## Chain 2:                0.780642 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'lm' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 3: Iteration: 2501 / 10000 [ 25%]  (Sampling)
## Chain 3: Iteration: 3500 / 10000 [ 35%]  (Sampling)
## Chain 3: Iteration: 4500 / 10000 [ 45%]  (Sampling)
## Chain 3: Iteration: 5500 / 10000 [ 55%]  (Sampling)
## Chain 3: Iteration: 6500 / 10000 [ 65%]  (Sampling)
## Chain 3: Iteration: 7500 / 10000 [ 75%]  (Sampling)
## Chain 3: Iteration: 8500 / 10000 [ 85%]  (Sampling)
## Chain 3: Iteration: 9500 / 10000 [ 95%]  (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.158748 seconds (Warm-up)
## Chain 3:                0.578621 seconds (Sampling)
## Chain 3:                0.737369 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'lm' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 4: Iteration: 2501 / 10000 [ 25%]  (Sampling)
## Chain 4: Iteration: 3500 / 10000 [ 35%]  (Sampling)
## Chain 4: Iteration: 4500 / 10000 [ 45%]  (Sampling)
## Chain 4: Iteration: 5500 / 10000 [ 55%]  (Sampling)
## Chain 4: Iteration: 6500 / 10000 [ 65%]  (Sampling)
## Chain 4: Iteration: 7500 / 10000 [ 75%]  (Sampling)
## Chain 4: Iteration: 8500 / 10000 [ 85%]  (Sampling)
## Chain 4: Iteration: 9500 / 10000 [ 95%]  (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.159318 seconds (Warm-up)
## Chain 4:                0.513896 seconds (Sampling)
## Chain 4:                0.673214 seconds (Total)
## Chain 4:
print(fit,pars=c('intercept','beta','sigma'),digits_summary = 3,
      probs = c(0.025,0.5,0.975))
## Inference for Stan model: lm.
## 4 chains, each with iter=10000; warmup=2500; thin=10; 
## post-warmup draws per chain=750, total post-warmup draws=3000.
## 
##            mean se_mean    sd   2.5%    50% 97.5% n_eff  Rhat
## intercept 0.000   0.001 0.035 -0.071 -0.001 0.067  2443 1.001
## beta      0.755   0.001 0.035  0.687  0.756 0.822  2922 1.000
## sigma     0.659   0.000 0.025  0.613  0.658 0.710  3125 1.001
## 
## Samples were drawn using NUTS(diag_e) at Thu Jun  3 17:52:56 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
       mean se_mean    sd   2.5%   50% 97.5% n_eff  Rhat

intercept 0.000 0.001 0.036 -0.071 0.001 0.073 2916 0.999 beta 0.755 0.001 0.035 0.689 0.755 0.824 2791 1.000 sigma 0.659 0.000 0.025 0.612 0.658 0.710 2974 0.999

compared to R’s lm above: intercept = 5.794e-17 (~ 0) weight = 7.547e-01 sigma = 0.6569515

plot(fit)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)