library(rjags)
library(tidyverse)
library(rjags)
library(ggmcmc)
library(polspline)
library(propagate)
intercept_mod <- "model{
# Priors
mu_int~dnorm(0, 0.0001) # Mean hyperparameter for random intercepts
sigma_int~dunif(0, 100) # SD hyperparameter for random intercepts
tau_int <- 1/(sigma_int*sigma_int)
for (i in 1:n_inds) {
alpha[i]~dnorm(mu_int, tau_int) # Random intercepts
}
beta~dnorm(0, 0.0001) # Common slope
sigma_res~dunif(0, 100) # Residual standard deviation
tau_res <- 1/(sigma_res*sigma_res)
# Likelihood
for (i in 1:n_obs) {
mu[i] <- alpha[inds[i]]+beta*x[i] # Expectation
y[i]~dnorm(mu[i], tau_res) # The actual (random) responses
}
}"
slope_model<- "model {
# Priors
mu_int~dnorm(0, 0.001) # Mean hyperparameter for random intercepts
sigma_int~dunif(0, 100) # SD hyperparameter for random intercepts
tau_int <- 1/(sigma_int*sigma_int)
mu_slope~dnorm(0, 0.001) # Mean hyperparameter for random slopes
sigma_slope~dunif(0, 100) # SD hyperparameter for slopes
tau_slope <- 1/(sigma_slope*sigma_slope)
for (i in 1:n_inds) {
alpha[i]~dnorm(mu_int, tau_int) # Random intercepts
beta[i]~dnorm(mu_slope, tau_slope) # Random slopes
}
sigma_res~dunif(0, 100) # Residual standard deviation
tau_res <- 1/(sigma_res*sigma_res) # Residual precision
# Likelihood
for (i in 1:n_obs) {
mu[i] <- alpha[inds[i]]+beta[inds[i]]*x[i]
y[i]~dnorm(mu[i], tau_res)
}
}"
Whales<-read.csv("https://tinyurl.com/aqm-data/whaleteeth.csv")
Whales$inds<-as.numeric(Whales$Whale)
dat<-list(
n_inds=length(unique(Whales$inds)),
n_obs =length(Whales$inds),
inds= Whales$inds,
y=Whales$X15N,
x=Whales$Age
)
model <- jags.model(textConnection(slope_model),
data = dat,
n.chains=1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 307
## Unobserved stochastic nodes: 27
## Total graph size: 1586
##
## Initializing model
update(model, 10000)
mod_sim <- coda.samples(model = model, variable.names = c("mu_int","mu_slope","alpha","beta"), n.iter = 10000)
summary(mod_sim)
##
## Iterations = 11001:21000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha[1] 13.16689 0.212665 2.127e-03 0.0060422
## alpha[2] 13.48769 0.300742 3.007e-03 0.0092333
## alpha[3] 12.50628 0.281675 2.817e-03 0.0073106
## alpha[4] 10.37204 0.255821 2.558e-03 0.0069144
## alpha[5] 12.47127 0.246334 2.463e-03 0.0065616
## alpha[6] 11.26491 0.270210 2.702e-03 0.0075070
## alpha[7] 11.06077 0.202957 2.030e-03 0.0052805
## alpha[8] 12.29949 0.336113 3.361e-03 0.0110099
## alpha[9] 13.29532 0.215434 2.154e-03 0.0055798
## alpha[10] 14.79853 0.274231 2.742e-03 0.0078438
## alpha[11] 11.77152 0.203833 2.038e-03 0.0057867
## beta[1] 0.05009 0.008660 8.660e-05 0.0002441
## beta[2] -0.01687 0.020192 2.019e-04 0.0006275
## beta[3] 0.10576 0.024527 2.453e-04 0.0006475
## beta[4] 0.16964 0.017831 1.783e-04 0.0004802
## beta[5] 0.04967 0.015947 1.595e-04 0.0004215
## beta[6] 0.10099 0.020602 2.060e-04 0.0005519
## beta[7] 0.15247 0.009342 9.342e-05 0.0002442
## beta[8] 0.10877 0.022907 2.291e-04 0.0007511
## beta[9] 0.10193 0.011043 1.104e-04 0.0002862
## beta[10] -0.06449 0.019967 1.997e-04 0.0005584
## beta[11] 0.11299 0.007723 7.723e-05 0.0002176
## mu_int 12.40457 0.474169 4.742e-03 0.0052126
## mu_slope 0.07903 0.026480 2.648e-04 0.0002921
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha[1] 12.74918 13.02198 13.16710 13.313376 13.57503
## alpha[2] 12.89247 13.27886 13.48496 13.689356 14.08042
## alpha[3] 11.94127 12.32216 12.50836 12.693707 13.04972
## alpha[4] 9.86746 10.20337 10.37658 10.540902 10.87770
## alpha[5] 11.98693 12.30602 12.47294 12.639464 12.94755
## alpha[6] 10.71313 11.08708 11.26906 11.447040 11.78700
## alpha[7] 10.66297 10.92483 11.06154 11.197103 11.46116
## alpha[8] 11.65636 12.06947 12.29477 12.524577 12.96059
## alpha[9] 12.87866 13.14859 13.29390 13.437187 13.72666
## alpha[10] 14.25869 14.62086 14.79873 14.979038 15.32925
## alpha[11] 11.36900 11.63393 11.77092 11.907895 12.17073
## beta[1] 0.03318 0.04415 0.05006 0.056006 0.06693
## beta[2] -0.05659 -0.03034 -0.01677 -0.003457 0.02248
## beta[3] 0.05799 0.08959 0.10565 0.122284 0.15481
## beta[4] 0.13508 0.15771 0.16936 0.181566 0.20545
## beta[5] 0.01880 0.03880 0.04975 0.060385 0.08098
## beta[6] 0.06076 0.08723 0.10070 0.114528 0.14219
## beta[7] 0.13402 0.14619 0.15240 0.158790 0.17079
## beta[8] 0.06321 0.09352 0.10905 0.124284 0.15278
## beta[9] 0.07977 0.09462 0.10203 0.109421 0.12345
## beta[10] -0.10416 -0.07766 -0.06444 -0.051534 -0.02394
## beta[11] 0.09773 0.10779 0.11308 0.118112 0.12808
## mu_int 11.44691 12.12023 12.41078 12.697286 13.35350
## mu_slope 0.02601 0.06242 0.07920 0.095327 0.13174
Musssels
d<-read.csv("https://tinyurl.com/aqm-data/mussels.csv")
str(d)
## 'data.frame': 113 obs. of 3 variables:
## $ Lshell : num 122.1 100.1 100.7 102.3 94.9 ...
## $ BTVolume: int 39 21 23 22 20 22 21 18 21 15 ...
## $ Site : Factor w/ 6 levels "Site_1","Site_2",..: 6 6 6 6 6 6 6 6 6 6 ...
summary(lm(data=d,Lshell~Site))
##
## Call:
## lm(formula = Lshell ~ Site, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.906 -8.340 1.031 9.231 30.550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.769 2.624 38.405 < 2e-16 ***
## SiteSite_2 8.467 3.748 2.259 0.0259 *
## SiteSite_3 6.037 5.409 1.116 0.2669
## SiteSite_4 -3.619 5.409 -0.669 0.5049
## SiteSite_5 18.697 3.925 4.763 6.02e-06 ***
## SiteSite_6 2.471 3.748 0.659 0.5111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.38 on 107 degrees of freedom
## Multiple R-squared: 0.2239, Adjusted R-squared: 0.1876
## F-statistic: 6.173 on 5 and 107 DF, p-value: 4.579e-05
data=list(y=d$Lshell,ind=as.numeric(d$Site),N=length(d$Lshell),p=length(levels(d$Site)))
modelstring="
model {
for (i in 1:N) {
mu[i]<-beta[ind[i]]
y[i] ~ dnorm(mu[i],tau)
}
for (j in 1:p) {
beta[j]~dnorm(0,0.0001)
for (n in 1:(j-1)){
difbeta[n,j]<-beta[n]-beta[j]
}
}
tau~dgamma(1,0.0001)
}
"
model=jags.model(textConnection(modelstring),data=data)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 113
## Unobserved stochastic nodes: 7
## Total graph size: 264
##
## Initializing model
update(model,n.iter=1000)
output=coda.samples(model=model,variable.names=c("difbeta"),n.iter=100000,thin=10)
summary(output)
##
## Iterations = 1010:101000
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## difbeta[1,2] -8.461 3.748 0.03748 0.03748
## difbeta[1,3] -5.826 5.429 0.05429 0.05820
## difbeta[2,3] 2.635 5.464 0.05464 0.05563
## difbeta[1,4] 3.721 5.362 0.05362 0.05362
## difbeta[2,4] 12.182 5.445 0.05445 0.05445
## difbeta[3,4] 9.547 6.711 0.06711 0.07059
## difbeta[1,5] -18.622 3.911 0.03911 0.03911
## difbeta[2,5] -10.161 3.951 0.03951 0.03951
## difbeta[3,5] -12.796 5.591 0.05591 0.05967
## difbeta[4,5] -22.343 5.561 0.05561 0.05561
## difbeta[1,6] -2.422 3.735 0.03735 0.03941
## difbeta[2,6] 6.039 3.811 0.03811 0.03811
## difbeta[3,6] 3.404 5.475 0.05475 0.05807
## difbeta[4,6] -6.143 5.430 0.05430 0.05430
## difbeta[5,6] 16.200 3.947 0.03947 0.03947
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## difbeta[1,2] -15.824 -10.9340 -8.450 -5.9865 -1.175
## difbeta[1,3] -16.508 -9.4872 -5.789 -2.2011 4.923
## difbeta[2,3] -8.115 -1.0113 2.703 6.2833 13.444
## difbeta[1,4] -6.733 0.1261 3.761 7.3489 14.272
## difbeta[2,4] 1.577 8.5929 12.192 15.7446 22.976
## difbeta[3,4] -3.601 5.0998 9.429 14.0164 22.798
## difbeta[1,5] -26.326 -21.2679 -18.629 -15.9639 -10.896
## difbeta[2,5] -17.835 -12.9231 -10.145 -7.4886 -2.394
## difbeta[3,5] -23.915 -16.4971 -12.784 -9.0562 -1.871
## difbeta[4,5] -33.591 -26.1326 -22.309 -18.6146 -11.245
## difbeta[1,6] -9.696 -4.9219 -2.407 0.1121 4.900
## difbeta[2,6] -1.410 3.5171 6.043 8.5554 13.579
## difbeta[3,6] -7.307 -0.2519 3.379 6.9903 14.406
## difbeta[4,6] -16.759 -9.7787 -6.146 -2.4707 4.453
## difbeta[5,6] 8.303 13.5803 16.216 18.8193 24.010
ms <-ggs(output)
mt<-filter(ms,grepl("difbeta",Parameter))
ggs_caterpillar(mt)
TukeyHSD(aov(d$Lshell~d$Site))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = d$Lshell ~ d$Site)
##
## $`d$Site`
## diff lwr upr p adj
## Site_2-Site_1 8.466769 -2.408905 19.342443 0.2201442
## Site_3-Site_1 6.037019 -9.660664 21.734702 0.8737518
## Site_4-Site_1 -3.619231 -19.316914 12.078452 0.9849444
## Site_5-Site_1 18.697436 7.305950 30.088922 0.0000867
## Site_6-Site_1 2.470769 -8.404905 13.346443 0.9859123
## Site_3-Site_2 -2.429750 -18.201132 13.341632 0.9976925
## Site_4-Site_2 -12.086000 -27.857382 3.685382 0.2355928
## Site_5-Site_2 10.230667 -1.262165 21.723498 0.1103764
## Site_6-Site_2 -5.996000 -16.977781 4.985781 0.6105029
## Site_4-Site_3 -9.656250 -29.069479 9.756979 0.7004668
## Site_5-Site_3 12.660417 -3.470986 28.791819 0.2123990
## Site_6-Site_3 -3.566250 -19.337632 12.205132 0.9862071
## Site_5-Site_4 22.316667 6.185264 38.448069 0.0015143
## Site_6-Site_4 6.090000 -9.681382 21.861382 0.8718474
## Site_6-Site_5 -16.226667 -27.719498 -4.733835 0.0011239