# R-Code zum Artikel 
# Kurze Einführung in Bayes-Statistik                                 
# Vogelwarte 3/2016                                                   

#-------------------------------------------------------------------------------------

# Abbildung 1                                                          
x <- 1:100
plot(x, dnorm(x, 50, 10)+0.004, type="l", lwd=2, axes=FALSE,
     xlab ="", ylab="Dichte", ylim=c(0, 0.05))
mtext("Parameterwert", side=1, line=3.5)
box(bty="L")
qn <- qnorm(p=c(0.01, 0.05, 0.95, 0.99), mean=50, sd=10)
weite <- 10
segments(c(0, qn[4]), c(0,0), c(qn[1], 100), c(0,0), col="red", lwd=weite, 
         lend="butt" )
segments(c(qn[1], qn[3]), c(0,0), c(qn[2], qn[4]), c(0,0), 
         col="violet", lwd=weite, lend="butt" )
segments(qn[2],0 , qn[3], 
         col="lightblue", lwd=weite, lend="butt" )
text(x=(c(0, qn)+c(qn,100))/2,y=rep(-0.005, 5), xpd=NA,
     col=c("red", "violet", "blue", "violet", "red"),
     labels=c("sehr\nwenig\nwahrscheinlich", "wenig\nwahrscheinlich", "wahrscheinlich", 
            "wenig\nwahrscheinlich", "sehr\nwenig\nwahrscheinlich"), 
     adj=c(0.5, 1))


# Beispiel Amsel
theta <- seq(0, 1, length=100)
priora <- 3
priorb <- 3
prior <- dbeta(theta, priora,priorb)
nmales <- 14
nfemales <- 10
qbeta(0.025, priora,priorb)
qbeta(0.975, priora,priorb)

likelihood <- dbinom(prob=theta, x=nmales, size=nmales+nfemales)/(sum(dbinom(prob=theta, x=nmales, size=nmales+nfemales))*theta[2])

posterior <- dbeta(theta, nmales+priora, nfemales+priorb)
# Mittelwert
(nmales+priora)/(nmales+priora+ nfemales+priorb)
qbeta(0.025, nmales+priora, nfemales+priorb)
qbeta(0.975, nmales+priora, nfemales+priorb)


# Abbildung 2                                                          
plot(0:1, range(posterior), type="n", xlab="Anteil Männchen", ylab="Dichte", las=1, cex.lab=1.4, cex.axis=1.2)
lines(theta, prior, lty=3, lwd=3)
lines(theta, likelihood, lty=2, lwd=3)
lines(theta, posterior, lty=1, lwd=3)
legend(0, 4, lty=c(3,2,1), lwd=3, legend=c("A-priori", "Daten", "A-posteriori"), bty="n")


#--------------------------------------------------------------------------------
# Beispiel Schneesperling                                              

geschlecht <- c("M", "M", "M", "W", "W", "W", "W", "W", "W")
fluegel <- c(125, 122, 121, 119, 119.5, 116, 121, 121, 117.5)
dat <- data.frame(geschlecht=geschlecht, fluegel=fluegel)
dat

mod <- lm(fluegel~geschlecht, data=dat)

                                                                
# Abbildung 3                                                  
par(mfrow=c(2,2))
plot(mod)

library(arm)
nsim <- 5000
bsim <- sim(mod, n.sim=nsim)
str(bsim)


# Abbildung 4                                                 
par(mfrow=c(1,2))
hist(bsim@coef[,1], main=names(coef(mod))[1], freq=FALSE, nclass=20)
hist(bsim@coef[,2], main=names(coef(mod))[2], freq=FALSE, nclass=20)

apply(bsim@coef, 2, mean)
apply(bsim@coef, 2, sd)
apply(bsim@coef, 2, quantile, prob=c(0.025, 0.975))

mw <- bsim@coef[,1] + bsim@coef[,2]
quantile(mw, prob=c(0.025, 0.975))


# Information aus dem Handbuch
mm <- 123.8
semm <- (128-121)/4/sqrt(21)


mw <- 116.7
semw <- (120-114)/4/sqrt(15)


mod <- lm(fluegel~geschlecht-1, data=dat)
mod

library(rstanarm)
stanmod <- stan_glm(fluegel~geschlecht-1, data=dat,
        prior = normal(c(mm, mw), scale=c(semm, semw)), 
        prior_intercept = NULL, iter=5000)

summary(stanmod)

bsim <- as.data.frame(stanmod)
d <-bsim$geschlechtW- bsim$geschlechtM
mean(d)

sum(d < -1)/nrow(bsim)

apply(bsim, 2, sd)

par(mfrow=c(1,3))
hist(bsim$geschlechtM)
hist(bsim$geschlechtW)
hist(bsim$sigma)



# Abbildung 5                                               
dat <- data.frame(Fall=LETTERS[1:5],
                    mean=c(1.55,1.55,0.25,0.25,0.25),
                    sd=c(0.25, 1, 1, 0.25, 0.1))
dat$lwr <- qnorm(p=0.025, mean=dat$mean, sd=dat$sd)
dat$upr <- qnorm(p=0.957, mean=dat$mean, sd=dat$sd)
dat$nht <- c("p<0.05", "ns", "ns", "ns", "p<0.05")
dat$postP <- round(1-pnorm(q=1, mean=dat$mean, sd=dat$sd) + pnorm(q= -1, mean=dat$mean, sd=dat$sd),2)

plot(1:5, seq(-1.8, 3.9, length=5), type="n", las=1, ylab="Effektgrösse",
     xlab=NA, xaxt="n", xlim=c(-0.5, 5.3), yaxt="n", mgp=c(2,0.6,0),
     cex.lab=1.2)
rect(-1, 1, 6, 5, border=NA, col="orange")
rect(-1, -6, 6, -1, border=NA, col="orange")
rect(-1, -1, 6, 1, border=NA, col="lightblue")
axis(2, tck=0.02, las=1, mgp=c(2, 0.6, 0))

abline(h=0)
segments(1:5, dat$lwr, 1:5, dat$upr, lwd=3, lend="butt")
points(1:5, dat$mean, pch=21, bg="white")
#text(1:5, dat$upr+0.25, dat$Fall)
axis(1, at=1:5, labels=dat$Fall, tck=-0.02)
text(c(1:5)-0.1, rep(3.8, 5), dat$nht, adj=c(0, NA))
text(-0.55, 3.8, "Nullhypothesetest", adj=c(0, NA))
text(c(1:5)-0.1, rep(3.45, 5), dat$postP, adj=c(0, NA))
text(-0.55, 3.45, "W. für Hypothese", adj=c(0, NA))

box()