## generating a dataframe on smoking, obesity, snoring and hypertonia ##

no.yes <- c("No","Yes")
smoking <- gl(2,1,8,no.yes)
obesity <- gl(2,2,8,no.yes)
snoring <- gl(2,4,8,no.yes)
n.tot <- c(60,17,8,2,187,85,51,23)
n.hyp <- c(5,2,1,0,35,13,15,8)
smoke.tab = data.frame(smoking,obesity,snoring,n.tot,n.hyp)

## n.hyp = number of patients with hypertonia, n.tot = all patients ##

## creating odds for the risk of hypertonia ##

hyp.tbl <- cbind(n.hyp,n.tot-n.hyp)

## important: 2nd column must fit the principle of odds ##

## not ##

logs = log(n.hyp/n.tot-n.hyp)

## because if all patients have hypertonia, we would divide by 0

## logistic regression model ##

h = glm(hyp.tbl~smoking+obesity+snoring,family=binomial("logit"))

## or ##

h = glm(hyp.tbl~smoking+obesity+snoring,family=binomial)

## alternative: proportions instead of logit ##

prop.hyp <- n.hyp/n.tot

## note: the total number of observations is missing for proportions. Thus: ##

glm.hyp <- glm(prop.hyp~smoking+obesity+snoring,binomial,weights=n.tot)

####### GLMM ########

## devoice.RData
## contains voiced/unvoiced responses for 10 listeners based on the proportion of voicing in alveolar fricatives ###

load("devoice.RData")
library(lme4)

barplot(table(devoice$prop.voice,devoice$response),beside=T)

d.v = devoice
d.v$response = droplevels(d.v$response,exclude="uv")

barplot(table(d.v$prop.voice,d.v$response),beside=T,names.arg=levels(d.v$prop.voice))


h = glmer(response~prop.voice+(1+prop.voice|subj),data=devoice,family="binomial")

## coef(h) contains the intercept k and the slope m for each participant

## the first element of the list is saved into a data.frame ##

d.coef = coef(h)$subj
names(d.coef) = c("k","m")

curve(exp(mean(d.coef$m)*x+mean(d.coef$k))/(1+exp(mean(d.coef$m)*x+mean(d.coef$k))),xlim=c(0,100),ylim=c(0,1))

## calculating the inflection point

ip.mean = mean(-d.coef$k/d.coef$m)

abline(v=ip.mean,col=4,lwd=2)