#################################################################
# Uebung zu Schaetzen und Testen 1 im WiSe 15/16 
# R-Loesung zur Uebung 12, Aufgabe 1
# Autor: Ludwig Bothmann; Aenderungen: David Ruegamer
#################################################################

# Funktionen laden
source("Rcode_blatt12_func.R")

# Datensatz einlesen
sprache <- read.table("http://www.stat.uni-muenchen.de/service/datenarchiv/sprache/sprache.asc", 
			  header=T)
str(sprache)

# Der Datensatz basiert auf einem Wahrnehmungsexperiments, 
# bei dem entschieden werden sollte, ob beim Sprechen eines Wortes 
# gelaechelt wurde oder nicht. Dazu wurden zwoelf verschiedene Woerter 
# jeweils acht Mal von einem Sprecher gesprochen und verschiedene
# physikalische Messgroessen hierzu ermittelt. Diese Woerter wurden 
# insgesamt 50 Versuchspersonen vorgespielt, die entscheiden mussten, 
# ob sie ein Laecheln hoerten oder nicht. Zusaetzlich war bekannt, 
# ob der Sprecher nach eigener Aussage beim Sprechen gelaechelt hatte.

# nummer  Referenznummer
# wort	gesprochenes Wort
# sprecher	Sprecher laechelt (1: ja, 0: nein)
# segtyp	kennzeichnet das Wortsegment, fuer das die physikalischen Messungen vorgenommen wurden
# meanf0	Mittelwert der Grundfrequenz
# slopef0	Steigung der Regressionsgeraden durch die Einzelmesswerte der Grundfrequenz
# meanf1	Frequenz des ersten Formanten
# meanf2	Frequenz des zweiten Formanten
# meanf3	Frequenz des dritten Formanten

# Versuchspersonen erkennen gut, ob der Sprecher laechelt, oder nicht
boxplot(anteil~sprecher, data=sprache, ylab="Laechelt", 
		  xlab="Anteil vorhergesagt",horizontal=T)

##################
# Aufgabe a)
##################

# Inwiefern kann das Laecheln durch gewisse physikalische Messgroessen
#	erkannt werden?
mod1 <- glm(sprecher ~ 1 + meanf0 + slopef0 + meanf1 + meanf2 + meanf3, 
				data=sprache, family=binomial)
summary(mod1)

# => meanf0, meanf1 und meanf3 zeigen einen signifikanten Zusammenhang

# Zusaetzlich die Einschaetzung der Versuchspersonen
mod2 <- glm(sprecher ~ 1 + anteil + meanf0 + slopef0 + meanf1 + meanf2 + meanf3, 
				data=sprache, family=binomial)
summary(mod2)
# "Zu gut" => Daten sind perfekt getrennt

# Korrelation zwischen Laecheln und Anteil der Personen, die auf Laecheln
#	getippt haben
cor(sprache$sprecher, sprache$anteil)


##################
# Aufgabe d)
##################

# Response-Vektor und Designmatrix
y <- sprache$sprecher
X <- as.matrix(cbind(1,sprache[,6:10]))
# Dimensionen
p <- ncol(X)
n <- length(y)

####################
# Priori-Parameter
####################

beta.0 <- rep(0,p)
B.0 <- diag(rep(1000, p))

####################
# Startwert fuer beta: beta.0
####################

beta.start <- t(beta.0) # beta.kq # beta.ml

# Anzahl Zufallszahlen aus der Posteriori
D <- 1000

##########################
# Metropolis-Algorithmus
##########################

set.seed(42)

bayes.out <- logit.sample(D=D,
								  y=y,
								  X=X,
								  B.0=B.0,
								  beta.0=beta.0,
								  beta.start=beta.start,
								  B.prop=vcov(mod1), # Kov-Matrix der Vorschlagsdichte
                                    # = Kov-Matrix der ML-S (bspw.)
								  do.for=TRUE) 

# Ergebnisse des Metropolis-Algorithmus
beta.sample <- bayes.out$beta.sample # Zufallszahlen aus Post von beta

# Wie viele Werte wurden akzeptiert?
d <- bayes.out$d 
d

# Plot der Samplingpfade
par(mfrow=c(2,3))
for(i in 1:p) plot(0:D, beta.sample[,i], type="l", main=names(coef(mod1))[i])
# => Konvergieren schnell. Man sieht gut, dass teilweise oft hintereinander
#	die neuen Vorschlaege nicht akzeptiert wurden.

# Die ersten burnin Samples werden verworfen
burnin <- 100


################################################
# Zum Vergleich: Als Vorschlags-Kov-Matrix B.kq
################################################

set.seed(42)

bayes.out.kq <- logit.sample(D=D,
									  y=y,
									  X=X,
									  B.0=B.0,
									  beta.0=beta.0,
									  beta.start=beta.start,
									  B.prop=vcov(lm(sprecher ~ 1 + meanf0 + slopef0 + meanf1 + meanf2 + meanf3, 
									                 data=sprache)),
									  do.for=TRUE)

# Ergebnisse des Metropolis-Algorithmus
beta.sample.kq <- bayes.out.kq$beta.sample
d.kq <- bayes.out.kq$d

d.kq
# => Es wird offensichtlich oefter akzeptiert als oben. Wenn man sich die
#	Samplingpfade ansieht merkt man, dass das kein Hinweis darauf ist,
#	dass das Sampling geglueckt ist.


# Plot der Samplingpfade
par(mfrow=c(2,3))
for(i in 1:p) plot(0:D, beta.sample.kq[,i], type="l", main=names(coef(mod1))[i])
# => Konvergieren ueberhaupt nicht, Vorschlaege anscheinend sehr schlecht,
#		einziger Unterschied: Kov-Matrix der Vorschlagsdichte


##################
# Aufgabe e)
##################

# Bayesianische Punktschaetzer: Arithmetisches Mittel der Samples
beta.bayes <- apply(beta.sample[-(1:burnin),], MARGIN=2, mean)
beta.bayes

# Vergleich mit KQ- und ML-Schaetzer
schaetzer <- round(cbind(coef(mod1), beta.bayes),4)
colnames(schaetzer) <- c("ML", "Bayes")
schaetzer
# => Schaetzung von Likelihood- und Bayes-Inferenz in etwa gleich 
#		in Richtung und Staerke

# Signifikanzniveau:
alpha <- 0.05

# Bayesianische Intervallschaetzer: Quantile der Samples
ki.bayes <- apply(beta.sample[-(1:burnin),], MARGIN=2, quantile, prob=c(alpha/2, 1-alpha/2))
colnames(ki.bayes) <- names(coef(mod1))
ki.bayes

# => meanf0, meanf1 und meanf3 signifikant, da die Kredibilitaetsintervalle
#		die 0 nicht ueberdecken
# => Wie bei ML: summary(mod1)

##################
# Aufgabe f)
##################

# Mit MCMClogit aus MCMCpack:

library(MCMCpack)

# ?MCMClogit

mclog <- MCMClogit(sprecher~1+meanf0 + slopef0 + meanf1 + meanf2 + meanf3,
						 data=sprache, burnin=100, mcmc=10000, thin=10,
						 beta.start=beta.0, b0=beta.0, B0=solve(B.0))
str(mclog)
plot(mclog, smooth=TRUE)


summary(mclog)

# Punktschaetzer
schaetzer2 <- round(cbind(schaetzer, summary(mclog)$statistics[,1]),4)
colnames(schaetzer2) <- c(colnames(schaetzer),"MCMClogit")
schaetzer2

# => Punktschaetzer in etwa gleich

# Intervallschaetzer (HPD und gleichendig)
round(t(HPDinterval(mclog, prob=0.95)),4)
ki.mcmclogit <- round(apply(mclog[,], MARGIN=2, quantile, prob=c(alpha/2, 1-alpha/2)),4)
ki.mcmclogit
# => Wieder meanf0, meanf1 und meanf3 signifikant

######################################
# Funktionen des Coda-Paketes
######################################

library(lattice)
library(coda)

help(package=coda)

# Eine Auswahl
acfplot(mclog) # Autokorrelationsplot
densityplot(mclog) # Dichteschaetzung der einzelnen beta_i
effectiveSize(mclog) # Effektive Laenge der jeweiligen Kette
							#	(Mit Autokorrelation verrechnet)
autocorr(mclog) # Autokorrelationen
crosscorr.plot(mclog) 	# Kreuzkorrelationen zwischen den Ketten der 
								#	einzelnen Kovariablen
levelplot(mclog) # Wie crosscor
HPDinterval(mclog, prob=1-alpha) # HPD-Intervalle
summary(mclog)$quantiles # Gleichendige Intervalle
qqmath(mclog) # QQ-Plots 
start(mclog) # Erster gespeicherter Wert
thin(mclog) # Ausduennungsintervall
time(mclog) # Zeitpunkt / Index der gespeicherten Werte
traceplot(mclog) # Plot der Pfade
xyplot(mclog) # Alternativer Plot der Pfade

# Konvergenzdiagnostik

# erzeuge zweite Kette, damit Varianzen zwischen den beiden
# Ketten und innerhalb der Ketten fuer Gelman-Rubin-Brooks-Plot
# verwendet werden koennen
mclog2 <- MCMClogit(sprecher~1+meanf0 + slopef0 + meanf1 + meanf2 + meanf3,
                   data=sprache, burnin=100, mcmc=10000, thin=10, seed=1337,
                   beta.start=beta.0, b0=beta.0, B0=solve(B.0))

mclogList <- mcmc.list(mclog,mclog2)
gelman.plot(mclogList)

##################
# Aufgabe h)
##################


###############
# exp(beta_1)
###############

# Punktschaetzer
mean(exp(beta.sample[-(1:burnin),2]))
mean(exp(mclog[,2]))
# => Sind alle anderen Kovariablen gleich, aendert sich die Chance auf 
#		Laecheln multiplikativ um 1.15 wenn der Wert von meanf0 um eine
#		Einheit steigt

# Intervallschaetzer
quantile(exp(beta.sample[-(1:burnin),2]), c(alpha/2, 1-alpha/2))
quantile(exp(mclog[,2]), c(alpha/2, 1-alpha/2))
# => Effekt der Kovariablen signifikant, da das KI die 1 nicht enthaelt

###################
# beta_1 + beta_2
###################

# Punktschaetzer
mean(beta.sample[-(1:burnin),2] + beta.sample[-(1:burnin),3])
mean(mclog[,2] + mclog[,3])

# Intervallschaetzer
quantile(beta.sample[-(1:burnin),2] + beta.sample[-(1:burnin),3], c(alpha/2, 1-alpha/2))
quantile(mclog[,2] + mclog[,3], c(alpha/2, 1-alpha/2))

# => Ergebnisse per Hand und per MCMClogit nahezu identisch