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])
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
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
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
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
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)
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)