City of Trees & BRMS

Trees

Estimating tree cover

In this excercise, we will be estimating the proportion of city area covered by trees in: Boise, ID.

but how much of the area of Boise is actually covered in trees?

Our first step will be to decide on an informed prior (you only need to do this once!). You will decide as a group. Use the beta distribution as your prior, and decide together on estimates for shape1 and shape2 that give a probability density that reflects what you think the distribution of tree cover in Boise is. In the example below, the estimated hyperpriors (or parameters in a prior distribution) are 10 and 5:

#Our group chose 7 & 3 as priors

curve(dbeta(x, 7, 3), col="red", main="To estimate/inform the prior")

Next, set up a grid sample from the posterior. Below, we are setting up a grid from 0.01 to 0.99, with a length of 1000 points. For each point we are calculating the likelihood of the beta distribution. Note how the plot is the same shape as the one above, but now is composed of discrete points (and remember, you will need to put in your own values for the beta distribution, not 0.5)

p_grid<-seq(from=0.01,to=0.99,length.out=1000)

prior<-dbeta(p_grid,7,3)
priorplot <- plot(prior~p_grid, col="purple")

#### Data collection

#use the code below to get a latitude random point and a longitude random point

#lat
runif(1,min=43.571530,max=43.626726)

## [1] 43.58544

#long
runif(1,min=-116.296842,max=-116.196592)

## [1] -116.266

Enter your random point in Google Earth. Then, identify whether your point is a tree (value of 1) or not a tree (value of 0). Enter your observation into the following Google Sheet.

Calculate likelihood

With data, in hand, you will now calculate the likelihood (all the relative ways the data could occur). In the example below, there were 10 total points collected, and 2 points with trees. You will, of course, need to update with your own values from the Google Sheet. In our class data there were 11 trees and 17 no trees.

likelihood<-dbinom(11,size=28,prob=p_grid)
likelihoodplot <- plot(likelihood~p_grid, col="steelblue")

In the code above, we are calculating the likelihood for each point in our grid, given the data. Note how the shape of the likelihood differs from the prior.

Normalize

The next step is to calculate the denominator in Bayes theorem–the marginal likelihood, or the sum of all the ways the data could occur. Dividing by this sum, will give us our posterior!

unnormalized.posterior<-likelihood*prior

posterior<-unnormalized.posterior/sum(unnormalized.posterior)

par(mfrow=c(1,3))

plot(prior~p_grid,main="Prior", col="purple")
plot(likelihood~p_grid,main="Likelihood", col="steelblue")
plot(posterior~p_grid,main="Posterior", col="coral")

We are done with our first iteration of Bayesian updating! Now, if you collect a second round of data, what can you use as your new prior?

YOU CAN USE posterior in as your new prior! that is Bayesian updating.

Summarize posterior

Now that we have successfully completed Bayesian updating, we can now summarize the posterior. The median is equivalent to a point estimate. But we can go beyond point estimates here! What if we want to describe the 95% probability tree cover is between two values?

samples<-sample(p_grid,prob=posterior,size=1e4,replace=TRUE)

median(samples)

## [1] 0.472042

quantile(samples,c(0.025,0.5,0.975))

##     2.5%      50%    97.5% 
## 0.318028 0.472042 0.629980

BRMS

Grants

Read in the data and clean it up. Let’s make Awarded = 1 and Declined = 0.

grants<-read.delim("GrantApplicationList.csv.txt")

#get rid of the pendings and withdrawns in the status column

pending <- subset(grants, Status=="Pending")
withdrawn <- subset(grants, Status=="Withdrawn")
awarded <- subset(grants, Status=="Awarded")
declined <- subset(grants, Status=="Declined")

#put it back together in dataframe
grants <- rbind(awarded, declined)

#give awardeds 1 and declined ZEROS YOOOO
grants$binary_Status <- ifelse(grants$Status=="Awarded", 1, 0)

Now, let’s change the year to years prior to 2025 that these were submitted. This will make the later years (closer to 2025) have a smaller number. We will call the new column year.

Notes on the code. dates in r are always a pain to work with. The code below seems complex, but we are basically taking the day we heard back about awarded or declined and then we are turning it into a date, because the .txt or whatever is not actually treating those values as dates. Next, we are turning those into integer values so we can subtract from 2025 and get a whole number. NOTE the “%Y” is to pull out JUST THE YEARS (we don’t want the month or day for this calculation )

grants$year <- 2025-as.integer(format(as.Date(grants$Status.Date, format = "%m/%d/%Y"), "%Y"))

Plot it!

plot(grants$year, grants$binary_Status, ylab="Awarded = 1, Declined = 0", xlab="Year")

A useful trick with binary data is to use jitter to add a small amount of random noise to the 0 or 1 data, to make it easier to see

plot(jitter(grants$binary_Status,0.5)~grants$year, ylab="Awarded = 1, Declined = 0", xlab="Year")

Model

\(y_i ~ Binomial(n_i, p_i)\)

or for single trials: \(y_i ~ Bernoulli(p_i)\)

Here the:

binary_Status = outcome \(y\)

year = predictor that influences \(p\)

bernoulli = binomial with n = 1

grant_model=brm(binary_Status~year,family=bernoulli,data=grants)

## Compiling Stan program...

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.2 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Trying to compile a simple C file

## Running "C:/PROGRA~1/R/R-45~1.2/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.3.0'
## gcc  -I"C:/PROGRA~1/R/R-45~1.2/include" -DNDEBUG   -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/Rcpp/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/src/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include "C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp"  -std=c++1y    -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include"      -O2 -Wall -std=gnu2x  -mfpmath=sse -msse2 -mstackrealign   -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
##                  from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
##                  from C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.2/etc/x64/Makeconf:289: foo.o] Error 1

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.2 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.009 seconds (Warm-up)
## Chain 1:                0.009 seconds (Sampling)
## Chain 1:                0.018 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.009 seconds (Warm-up)
## Chain 2:                0.008 seconds (Sampling)
## Chain 2:                0.017 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.009 seconds (Warm-up)
## Chain 3:                0.008 seconds (Sampling)
## Chain 3:                0.017 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.008 seconds (Warm-up)
## Chain 4:                0.009 seconds (Sampling)
## Chain 4:                0.017 seconds (Total)
## Chain 4:

Note that there are other ways we could write this!

brm(binary_Status~year,family=binomial,data=grants)

brm(binary_Status~year,family=binomial(link=logit),data=grants)

Extract Posterior Samples

posterior_samples <- as.data.frame(grant_model)

Question about the probability of direction: what is the probability the effect is positive or negative?

Remember: Bayesian analysis is counting!

length(which(posterior_samples$b_year>0))/length(posterior_samples$b_year)

## [1] 0.90725

Plot it

histogram <- hist(posterior_samples$b_year,breaks=100, main="Probability of Direction", xlab="Posterior b_year")

histogram

## $breaks
##   [1] -0.39 -0.38 -0.37 -0.36 -0.35 -0.34 -0.33 -0.32 -0.31 -0.30 -0.29 -0.28
##  [13] -0.27 -0.26 -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16
##  [25] -0.15 -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04
##  [37] -0.03 -0.02 -0.01  0.00  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08
##  [49]  0.09  0.10  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.20
##  [61]  0.21  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.30  0.31  0.32
##  [73]  0.33  0.34  0.35  0.36  0.37  0.38  0.39  0.40  0.41  0.42  0.43  0.44
##  [85]  0.45  0.46  0.47  0.48  0.49  0.50  0.51  0.52  0.53  0.54  0.55  0.56
##  [97]  0.57  0.58  0.59  0.60  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68
## [109]  0.69  0.70  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.80
## [121]  0.81  0.82  0.83  0.84
## 
## $counts
##   [1]   1   0   0   1   1   1   0   0   0   0   1   0   0   0   2   2   4   2
##  [19]   1   6   1   9   5  14  10   9  13  11  12  12  22  21  21  20  18  28
##  [37]  45  38  40  49  46  43  67  64  61  55  56  70  81  76  79  86  94  90
##  [55] 100  85  79 101 105  82  94  97  94  76  85  89  80  91  81  92  79  71
##  [73]  76  62  63  54  59  72  46  48  59  37  29  39  29  29  23  30  26  28
##  [91]  19  16  15  14  23   9  15   9  12   8   3   9   6   7  11   9   6   4
## [109]   2   2   2   6   4   4   1   0   1   0   0   2   0   1   2
## 
## $density
##   [1] 0.025 0.000 0.000 0.025 0.025 0.025 0.000 0.000 0.000 0.000 0.025 0.000
##  [13] 0.000 0.000 0.050 0.050 0.100 0.050 0.025 0.150 0.025 0.225 0.125 0.350
##  [25] 0.250 0.225 0.325 0.275 0.300 0.300 0.550 0.525 0.525 0.500 0.450 0.700
##  [37] 1.125 0.950 1.000 1.225 1.150 1.075 1.675 1.600 1.525 1.375 1.400 1.750
##  [49] 2.025 1.900 1.975 2.150 2.350 2.250 2.500 2.125 1.975 2.525 2.625 2.050
##  [61] 2.350 2.425 2.350 1.900 2.125 2.225 2.000 2.275 2.025 2.300 1.975 1.775
##  [73] 1.900 1.550 1.575 1.350 1.475 1.800 1.150 1.200 1.475 0.925 0.725 0.975
##  [85] 0.725 0.725 0.575 0.750 0.650 0.700 0.475 0.400 0.375 0.350 0.575 0.225
##  [97] 0.375 0.225 0.300 0.200 0.075 0.225 0.150 0.175 0.275 0.225 0.150 0.100
## [109] 0.050 0.050 0.050 0.150 0.100 0.100 0.025 0.000 0.025 0.000 0.000 0.050
## [121] 0.000 0.025 0.050
## 
## $mids
##   [1] -0.385 -0.375 -0.365 -0.355 -0.345 -0.335 -0.325 -0.315 -0.305 -0.295
##  [11] -0.285 -0.275 -0.265 -0.255 -0.245 -0.235 -0.225 -0.215 -0.205 -0.195
##  [21] -0.185 -0.175 -0.165 -0.155 -0.145 -0.135 -0.125 -0.115 -0.105 -0.095
##  [31] -0.085 -0.075 -0.065 -0.055 -0.045 -0.035 -0.025 -0.015 -0.005  0.005
##  [41]  0.015  0.025  0.035  0.045  0.055  0.065  0.075  0.085  0.095  0.105
##  [51]  0.115  0.125  0.135  0.145  0.155  0.165  0.175  0.185  0.195  0.205
##  [61]  0.215  0.225  0.235  0.245  0.255  0.265  0.275  0.285  0.295  0.305
##  [71]  0.315  0.325  0.335  0.345  0.355  0.365  0.375  0.385  0.395  0.405
##  [81]  0.415  0.425  0.435  0.445  0.455  0.465  0.475  0.485  0.495  0.505
##  [91]  0.515  0.525  0.535  0.545  0.555  0.565  0.575  0.585  0.595  0.605
## [101]  0.615  0.625  0.635  0.645  0.655  0.665  0.675  0.685  0.695  0.705
## [111]  0.715  0.725  0.735  0.745  0.755  0.765  0.775  0.785  0.795  0.805
## [121]  0.815  0.825  0.835
## 
## $xname
## [1] "posterior_samples$b_year"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

col=ifelse(histogram$mids < 0, "steelblue", "tomato")

hist(posterior_samples$b_year,breaks=100,col=col, main="Probability of Direction", xlab="Posterior b_year")

This is a histogram that shows the probability of direction, graphically b_year tells you how much the log-odds of getting a grant changes for a 1-unit increase in year.

Now color it so that values below 0 are blue and others are tomatoes!

b_year = 0 \(\rightarrow\) no relationship between year and grant success

b_year > 0 \(\rightarrow\) positive effect

b_year < 0 \(\rightarrow\) negative effect

Extract Samples Directly!

Here we will create some vectors of data taking the slopes and year values out of the grant model.

slope_values<-posterior_samples$b_year

intercept_values<-posterior_samples$b_Intercept

Plot it!

hist(intercept_values,breaks=100)

We can look at the median, our single “best guess” at the parameter values.

median(slope_values)

## [1] 0.2170567

median(intercept_values)

## [1] -1.967874

This is the probability that the effect is between two numbers on a given interval. We can start by looking at the standard 95% CI

quantile(slope_values,c(0.025,0.975))

##       2.5%      97.5% 
## -0.1058066  0.5921165

But maybe 90% is better suited for this question, given that we know there is a lot of randomness with getting grants.

quantile(slope_values,c(0.1,0.9))

##         10%         90% 
## 0.006978962 0.446004491

Now let’s make a nice plot showing posterior effects:

plot(jitter(grants$binary_Status,0.5)~grants$year,pch="$",col="green4", main="Posterior Effects", xlab="Status", ylab="Year")

curve(plogis(median(intercept_values)+median(slope_values)*x),add=T,lwd=5)


for(i in 1:length(intercept_values)){
  curve(plogis(intercept_values[i]+slope_values[i]*x),add=T,lwd=0.5,col=rgb(0, 0, 1, alpha = 0.4))
}

curve(plogis(median(intercept_values)+median(slope_values)*x),add=T,lwd=5)

Above we made a big clunky for loop. This is Trevor’s preference, but many other ways to do this!

City of Trees MODEL

Here is another model in brms, that translates to our city of trees activity. Put in your own values for the data and prior, testing for trees:

y=rbinom(n=5,size=1,prob=0.9)

dat=data.frame(y=y)

#note that our link is an identity link. this means that the probability
#is not transformed (there is no logit link). How does this mirror what we did in class?


treemodel <-brm(y ~ 1,data=dat, family=bernoulli(link="identity"), prior  = prior(beta(1, 1), class = "Intercept",lb=0,ub=1))

## Compiling Stan program...

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.2 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Trying to compile a simple C file

## Running "C:/PROGRA~1/R/R-45~1.2/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.3.0'
## gcc  -I"C:/PROGRA~1/R/R-45~1.2/include" -DNDEBUG   -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/Rcpp/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/src/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/"  -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/jessi/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include "C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp"  -std=c++1y    -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include"      -O2 -Wall -std=gnu2x  -mfpmath=sse -msse2 -mstackrealign   -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
##                  from C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
##                  from C:/Users/jessi/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
##                  from <command-line>:
## C:/Users/jessi/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
##   679 | #include <cmath>
##       |          ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.2/etc/x64/Makeconf:289: foo.o] Error 1

## WARNING: Rtools is required to build R packages, but is not currently installed.
## 
## Please download and install the appropriate version of Rtools for 4.5.2 from
## https://cran.r-project.org/bin/windows/Rtools/.

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.007 seconds (Warm-up)
## Chain 1:                0.005 seconds (Sampling)
## Chain 1:                0.012 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 3e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.006 seconds (Warm-up)
## Chain 2:                0.006 seconds (Sampling)
## Chain 2:                0.012 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.005 seconds (Warm-up)
## Chain 3:                0.004 seconds (Sampling)
## Chain 3:                0.009 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 4e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.005 seconds (Warm-up)
## Chain 4:                0.005 seconds (Sampling)
## Chain 4:                0.01 seconds (Total)
## Chain 4:

#the prior code specifies a beta prior for the intercept. it works better

#when we specify a lower bound (lb) of 0 and an upperbound (ub) of 1, even though

#this is kind of already baked into the beta

summary(treemodel)

##  Family: bernoulli 
##   Links: mu = identity 
## Formula: y ~ 1 
##    Data: dat (Number of observations: 5) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.86      0.12     0.55     1.00 1.00     1006      908
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).