Chapter12
Easy
12E1
The prior with smaller standard deviation will produce the most shrinkage.
12E2
alist(y~dbinom(1,p),
logit(p) <- a[group] + b*X,
a[group] ~ dnorm(0,10),
b ~ dnorm(0,5))12E3
alist(y~dnorm(mu,sigma),
mu=a[group]+ b*x,
a[group] ~ dnorm(a, sigma_a),
b ~ dnorm(0,1),
sigma ~ dcauchy(0,1),
sigma_a ~ dcauchy(0,1))
# curve(dcauchy(x, 0, 1), from = 0, to = 10)12E4
alist(y~dpois(lambda),
log(lambda) = a + b[group]*x,
a ~ dnorm(0,10),
b[group] ~ dnorm(b,sigma_b),
sigma_b ~ dcauchy(0,1))12E5
alist(y~dpois(lambda),
log(lambda) = a +a_group[group] +a_sex[sex] + b*x,
a_group[group] ~ dnorm(0,sigma_a_group),
a_sex[sex] ~ dnorm(0,sigma_a_sex),
b ~ dnorm(0,5),
c(sigma_a_group,sigma_a_sex) ~ dcauchy(0,1) )Medium
12M1
data(reedfrogs)
# dummy variables
reedfrogs$pred <- ifelse(reedfrogs$pred=="pred",1,0)
reedfrogs$frogsize <- ifelse(reedfrogs$size=="big",1,0)
reedfrogs$tank <- 1:nrow(reedfrogs)
# m12.2 from the book
mm_i <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f413b4339e.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.229 seconds (Warm-up)
## 0.8 seconds (Sampling)
## 2.029 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.044 seconds (Warm-up)
## 1.045 seconds (Sampling)
## 2.089 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.091 seconds (Warm-up)
## 0.749 seconds (Sampling)
## 1.84 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.296 seconds (Warm-up)
## 2.072 seconds (Sampling)
## 3.368 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.016 seconds (Sampling)
## 0.016 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]mm_ip <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] + bP*pred ,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1),
bP ~ dnorm(0,5)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f4359c63d6.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.53 seconds (Warm-up)
## 1.314 seconds (Sampling)
## 2.844 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.375 seconds (Warm-up)
## 1.357 seconds (Sampling)
## 2.732 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.477 seconds (Warm-up)
## 1.233 seconds (Sampling)
## 2.71 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.45 seconds (Warm-up)
## 1.309 seconds (Sampling)
## 2.759 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]mm_is <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] +bS*frogsize,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1),
bS ~ dnorm(0,5)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f475585a0e.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 4.385 seconds (Warm-up)
## 1.689 seconds (Sampling)
## 6.074 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.64 seconds (Warm-up)
## 1.44 seconds (Sampling)
## 3.08 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.531 seconds (Warm-up)
## 1.369 seconds (Sampling)
## 2.9 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.695 seconds (Warm-up)
## 1.522 seconds (Sampling)
## 3.217 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]mm_ips <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] + bP*pred + bS*frogsize,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1),
c(bS,bP) ~ dnorm(0,5)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f45ee917df.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.597 seconds (Warm-up)
## 1.377 seconds (Sampling)
## 2.974 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.487 seconds (Warm-up)
## 1.453 seconds (Sampling)
## 2.94 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.665 seconds (Warm-up)
## 1.414 seconds (Sampling)
## 3.079 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 1.583 seconds (Warm-up)
## 1.441 seconds (Sampling)
## 3.024 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]mm_ipsi <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] + bP*pred + bS*frogsize + bPS*pred*frogsize,
a_tank[tank] ~ dnorm( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1),
c(bS,bP,bPS) ~ dnorm(0,5)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f44191780.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 2.222 seconds (Warm-up)
## 2.593 seconds (Sampling)
## 4.815 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 2.14 seconds (Warm-up)
## 2.619 seconds (Sampling)
## 4.759 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 2.494 seconds (Warm-up)
## 2.718 seconds (Sampling)
## 5.212 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 2.309 seconds (Warm-up)
## 2.177 seconds (Sampling)
## 4.486 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]compare(mm_i,mm_ip,mm_is,mm_ips, mm_ipsi)## WAIC pWAIC dWAIC weight SE dSE
## mm_ip 1001.4 28.8 0.0 0.35 37.31 NA
## mm_ips 1001.5 28.0 0.1 0.33 37.34 1.69
## mm_ipsi 1001.6 27.9 0.2 0.31 37.56 2.99
## mm_i 1010.3 38.0 8.9 0.00 37.94 6.71
## mm_is 1010.4 38.1 9.0 0.00 38.03 6.79all models containing “predator” are equally good
pr_ip <- precis(mm_ip , depth=2)
pr_ips <- precis(mm_ips , depth=2)
pr_ipsi <- precis(mm_ipsi, depth=2)
plot(y=reedfrogs$propsurv, x=reedfrogs$tank, ylim=c(0,1), pch=19)
abline(h=median(reedfrogs$propsurv))
points(y=inv_logit(pr_ip@output$Mean[1:48]) , x=reedfrogs$tank, col="red")
points(y=inv_logit(pr_ips@output$Mean[1:48]) , x=reedfrogs$tank, col="green")
points(y=inv_logit(pr_ipsi@output$Mean[1:48]), x=reedfrogs$tank, col="blue")
Adding the predator predictor raises all the intercepts to above the mean, probably the presence of predators is the thing that causes poor survival in some tanks.
12M2
post_ip <- extract.samples(mm_ip)
post_ips <- extract.samples(mm_ips)
post_ipsi <- extract.samples(mm_ipsi)
post_is <- extract.samples(mm_is)
post_i <- extract.samples(mm_i)
par(mfrow=c(2,2))
dens(post_ip$bP, main="b_predator",xlim=c(-5,1))
dens(post_ips$bP, add=T, col="red")
dens(post_ipsi$bP, add=T, col="blue")
dens(post_is$bS, main="b_size",xlim=c(-3,3), ylim=c(0,1.2))
abline(v=0)
dens(post_ips$bS, add=T, col="red")
dens(post_ipsi$bS, add=T, col="blue")
dens(post_ipsi$bPS, main="b_interaction", col="blue")
plot.new()
legend("topleft",legend=c("single variable","both variables","interaction"),fill=c("black","red","blue"))
Predator has a strong negative effect on survival, while the effect of size is modest. All of the models containing predator fit equally well. For the interaction model, the effect of predator seems divided over the single and interaction coefficients.
12M3
mm_i_cauchy <- map2stan(
alist(
surv ~ dbinom( density , p ) ,
logit(p) <- a_tank[tank] ,
a_tank[tank] ~ dcauchy( a , sigma ) ,
a ~ dnorm(0,1) ,
sigma ~ dcauchy(0,1)
), data=reedfrogs , iter=4000 , chains=4 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f4e5a3cdd.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 4.162 seconds (Warm-up)
## 10.989 seconds (Sampling)
## 15.151 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 26.315 seconds (Warm-up)
## 30.866 seconds (Sampling)
## 57.181 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 4.357 seconds (Warm-up)
## 21.923 seconds (Sampling)
## 26.28 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 4000 [ 0%] (Warmup)
## Iteration: 400 / 4000 [ 10%] (Warmup)
## Iteration: 800 / 4000 [ 20%] (Warmup)
## Iteration: 1200 / 4000 [ 30%] (Warmup)
## Iteration: 1600 / 4000 [ 40%] (Warmup)
## Iteration: 2000 / 4000 [ 50%] (Warmup)
## Iteration: 2001 / 4000 [ 50%] (Sampling)
## Iteration: 2400 / 4000 [ 60%] (Sampling)
## Iteration: 2800 / 4000 [ 70%] (Sampling)
## Iteration: 3200 / 4000 [ 80%] (Sampling)
## Iteration: 3600 / 4000 [ 90%] (Sampling)
## Iteration: 4000 / 4000 [100%] (Sampling)
##
## Elapsed Time: 7.72 seconds (Warm-up)
## 32.258 seconds (Sampling)
## 39.978 seconds (Total)
##
##
## SAMPLING FOR MODEL 'surv ~ dbinom(density, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]pr_i <- precis(mm_i, depth = 2)
pr_i_cauchy <- precis(mm_i_cauchy, depth=2)
plot(y=reedfrogs$propsurv, x=reedfrogs$tank, ylim=c(0,1), pch=19)
abline(h=median(reedfrogs$propsurv))
points(y=inv_logit(pr_i@output$Mean[1:48]) , x=reedfrogs$tank, col="red")
points(y=inv_logit(pr_i_cauchy@output$Mean[1:48]) , x=reedfrogs$tank, col="blue")
legend("bottomleft", legend=c("original data", "a_sigma ~ dnorm(0,1) ","a_sigma ~ dcauchy(0,1)"),fill=c("black","red","blue"))
The cauchy prior does not always shrink the tank means toward the median. Does it overfit more because it allows extreme values for the intercepts?
12M4
rm(list=ls())
data("chimpanzees")
d <- chimpanzees
d$recipient <- NULL # get rid of NAs
d$block_id <- d$block # name 'block' is reserved by Stan
# original
m12.5 <- map2stan(
alist(
pulled_left ~ dbinom( 1 , p ),
logit(p) <- a + a_actor[actor] + a_block[block_id] +
(bp + bpc*condition)*prosoc_left,
a_actor[actor] ~ dnorm( 0 , sigma_actor ),
a_block[block_id] ~ dnorm( 0 , sigma_block ),
c(a,bp,bpc) ~ dnorm(0,10),
sigma_actor ~ dcauchy(0,1),
sigma_block ~ dcauchy(0,1)
) ,
data=d, warmup=1000 , iter=6000 , chains=4 , cores=3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f479692285.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'pulled_left ~ dbinom(1, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 2000 / 20000 ]
[ 4000 / 20000 ]
[ 6000 / 20000 ]
[ 8000 / 20000 ]
[ 10000 / 20000 ]
[ 12000 / 20000 ]
[ 14000 / 20000 ]
[ 16000 / 20000 ]
[ 18000 / 20000 ]
[ 20000 / 20000 ]# alternative
m12.5_alt <- map2stan(
alist(
pulled_left ~ dbinom( 1 , p ),
logit(p) <- a_actor[actor] + a_block[block_id] +
(bp + bpc*condition)*prosoc_left,
a_actor[actor] ~ dnorm( a , sigma_actor ),
a_block[block_id] ~ dnorm( g , sigma_block ),
c(a,g,bp,bpc) ~ dnorm(0,10),
sigma_actor ~ dcauchy(0,1),
sigma_block ~ dcauchy(0,1)
) ,
data=d, warmup=1000 , iter=6000 , chains=4 , cores=3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file24f447b75e62.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'pulled_left ~ dbinom(1, p)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 2000 / 20000 ]
[ 4000 / 20000 ]
[ 6000 / 20000 ]
[ 8000 / 20000 ]
[ 10000 / 20000 ]
[ 12000 / 20000 ]
[ 14000 / 20000 ]
[ 16000 / 20000 ]
[ 18000 / 20000 ]
[ 20000 / 20000 ]precis(m12.5, depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_actor[1] -1.15 0.99 -2.58 0.35 3119 1
## a_actor[2] 4.17 1.58 1.92 6.43 5814 1
## a_actor[3] -1.46 0.99 -2.91 0.01 3092 1
## a_actor[4] -1.45 0.99 -2.87 0.05 3116 1
## a_actor[5] -1.15 0.99 -2.63 0.29 3107 1
## a_actor[6] -0.20 0.99 -1.73 1.21 3103 1
## a_actor[7] 1.34 1.02 -0.16 2.86 3322 1
## a_block[1] -0.18 0.23 -0.54 0.12 4817 1
## a_block[2] 0.04 0.19 -0.25 0.33 11394 1
## a_block[3] 0.05 0.19 -0.23 0.35 8737 1
## a_block[4] 0.01 0.18 -0.29 0.30 12190 1
## a_block[5] -0.03 0.19 -0.33 0.26 11980 1
## a_block[6] 0.11 0.21 -0.17 0.44 6425 1
## a 0.43 0.98 -1.03 1.86 3018 1
## bp 0.83 0.26 0.42 1.26 11394 1
## bpc -0.14 0.30 -0.61 0.34 11255 1
## sigma_actor 2.26 0.93 1.06 3.41 4487 1
## sigma_block 0.23 0.18 0.01 0.43 2320 1precis(m12.5_alt, depth=2)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_actor[1] 1.42 7.36 -9.96 9.52 5 1.29
## a_actor[2] 6.55 7.30 -5.28 14.39 5 1.26
## a_actor[3] 1.17 7.41 -10.15 9.39 5 1.30
## a_actor[4] 1.17 7.41 -10.18 9.26 5 1.30
## a_actor[5] 1.43 7.36 -9.94 9.51 5 1.29
## a_actor[6] 2.39 7.37 -8.93 10.47 5 1.29
## a_actor[7] 3.93 7.36 -7.54 11.96 5 1.29
## a_block[1] -2.32 7.31 -10.15 9.27 5 1.28
## a_block[2] -2.11 7.37 -10.00 9.44 5 1.30
## a_block[3] -2.10 7.38 -10.05 9.37 5 1.30
## a_block[4] -2.15 7.36 -10.18 9.26 5 1.29
## a_block[5] -2.18 7.35 -10.06 9.39 5 1.29
## a_block[6] -2.05 7.39 -9.96 9.47 5 1.30
## a 2.61 7.46 -8.42 11.23 5 1.30
## g -2.15 7.36 -10.00 9.45 5 1.29
## bp 0.83 0.23 0.51 1.25 34 1.06
## bpc -0.20 0.28 -0.73 0.18 45 1.07
## sigma_actor 2.39 0.81 1.22 3.48 54 1.05
## sigma_block 0.21 0.19 0.03 0.43 6 1.25In the first model, a_actor and a_block are divergence from the mean a, and independent of eachother. In the second model, this is no longer the case, because their means are not zero. This means that there are fewer effective samples to determine their values.
post_ori <- extract.samples(m12.5)
post_alt <- extract.samples(m12.5_alt)
par(mfrow=c(1,2))
for(i in 1:length(unique(d$actor))){
dens(post_ori$a_actor[,i], add=(i!=1), main="a_actor", xlim=c(-20,20) )
dens(post_alt$a_actor[,i], add=T, col="red")
}
for(i in 1:length(unique(d$block_id))){
dens(post_ori$a_block[,i], add=(i!=1), main="a_block", xlim=c(-20,20) )
dens(post_alt$a_block[,i], add=T, col="red")
}
legend("topright",legend=c("original","alternative"),fill=c("black","red"))
This is also reflected in the posterior distributions: all alphas become very flat and wide, because a_actor depends on a_block.
compare(m12.5,m12.5_alt)## WAIC pWAIC dWAIC weight SE dSE
## m12.5_alt 530.9 9.4 0.0 0.69 19.47 NA
## m12.5 532.6 10.4 1.6 0.31 19.67 0.68Still, the models have near identical WAIC values. However, the alternative model is near impossible to interpret, so it is still a bad idea to make models this way.
Hard
12H1
rm(list=ls())
data(bangladesh)
d <- bangladesh
d$district_id <- as.integer(as.factor(d$district))
data <- list(district_id = d$district_id, use_contraception = d$use.contraception)
# will not run as markdown
modfile="models12H1.Rdata"
if(file.exists(modfile)){load(modfile)}else{
# fixed effects model, fitting a different alpha for each district
m12H1.fe <- map2stan(
alist(
use.contraception ~ dbinom(1, p) ,
logit(p) <- alpha[district_id],
alpha[district_id] ~ dnorm(0, 10)
), data=d )
# multilevel model, alpha per district depends on overall distribution of alpha
m12H1.ml <- map2stan(
alist(
use.contraception ~ dbinom(1, p) ,
logit(p) <- alpha[district_id],
alpha[district_id] ~ dnorm(a, sigma),
a ~ dnorm(0, 10),
sigma ~ dcauchy(0, 1)
), data=d )
save(m12H1.fe,m12H1.ml, file =modfile)
}link_fe <- link(m12H1.fe, data=data.frame(district_id=1:60))## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]link_fe_mean <- apply(link_fe,2,mean)
link_fe_PI <- apply(link_fe,2,PI)
link_ml <- link(m12H1.ml, data=data.frame(district_id=1:60))## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]link_ml_mean <- apply(link_ml,2,mean)
link_ml_PI <- apply(link_ml,2,PI)
# get proportions of original
library(dplyr)
d_sum <- d %>% as.tbl %>% group_by(district_id) %>% summarize(p_use_contraception=mean(use.contraception))
plot(y=d_sum$p_use_contraception, x=d_sum$district_id, pch=19, xlab="districts", ylab="proportion using contraception")
abline(h=median(d_sum$p_use_contraception))
legend("topright",legend=c("data","fixed effects", "multilevel"),fill=c("black","red","blue"))
points(y=link_fe_mean,x=1:60, col="red")
for(i in 1:60){segments(x0=i,x1=i,y0=link_fe_PI[1,i],y1=link_fe_PI[2,i],col="red")}
points(y=link_ml_mean,x=(1:60)+0.2, col="blue")
for(i in 1:60){segments(x0=i+0.2,x1=i+0.2,y0=link_ml_PI[1,i],y1=link_ml_PI[2,i],col="blue")}
The fixed effects model closely follows the data in terms of means (the mean is not accurate when the true proportion is 0 or 1).
The mixed effects models shrinks all extreme points to the median.
12H2
Stan doesn’t like this a1,a2,a3,a4,a5,a6 vector, I can’t figure out how to make it work. Cannot fit the ML model until I figure this out.
12H3
Same problem