##                                Mandag 11. jan

x <- 3
log(x+1)
# Vektoriserte beregninger
vekt <- c(60, 72, 57, 90, 95, 72)
logvekt <- log(vekt)
høyde <- c(1.75, 1.80, 1.65, 1.9, 1.74, 1.91)
bmi <- vekt/høyde^2 # 2-tallet (en vektor av lengde 1 resirkuleres)

# egenskaper ved en vektor
sum(vekt)
n <- length(vekt)
var(vekt)
sd(vekt)

# standardavviket beregnet uten bruk av sd
sqrt(sum((vekt - mean(vekt))^2)/(length(vekt)-1))

# enkle statistiske tester
t.test(bmi, mu=22.5) # Ett-utvalgs t-test
minmodell <- lm(vekt ~ høyde) # Lineær regresjon
residuals(minmodell) # Videreberegninger på modellobjekte
coef(minmodell)
summary(minmodell)

# Grafikk
plot(høyde, vekt)
coef(minmodell)
abline(a=-46.33, b=67.35)
abline(minmodell)
points(høyde, 22.5*høyde^2,col="red",pch="+")
curve(sin(x),0,2*pi)

# Definering av egne funksjoner.  Funksjonen returnerer verdien av siste uttrykk i funksjonskroppen.
areal.av.sirkel <- function(r) {
  2*pi*r^2
}
areal.av.sirkel(10)
areal.av.sirkel(5)

# Vektorer
bmi > 25 # Den logiske operatoren ">" virker elementvis på bmi og resultatet blir en logisk vektor
seq(1,10,length=20) # 
seq(1,10,by=.5)
1:10 # Det samme som seq(1,10,by=1)
rep(5,10) 

# Matriser (og linear algebra for de som har hatt MA0002 eller BI2033)
matrix(1:12,3,4) # Elementene fylles inn søylevis som default
A <- matrix(c(0, 2, 3, 1, # En Leslie-matrise
              .5,0, 0, 0,
              0,.6,0,0,
              0,0,.2,0), nrow=4,ncol=4, byrow=TRUE)

#                                  Torsdag 14.jan. 

A <- rbind(c(0,2,3,1),cbind(diag(c(.5,.6,.2)),rep(0,3))) # Samme Lesliematrise (se ?cbind,rbind og diag)
eigen(A)$val[1] # Første egenverdi (populasjonens vekstrate)
vec <- eigen(A)$vec[,1] # Første egenvektor
vec/sum(vec) # er den stabile aldersfordelingen
N0 <- c(10,0,0,0) # Populasjonen ved t=0
A %*% N0 # Populasjonen ved t=1,2,3 beregnet ved matrisemultiplikasjon %*%
A %*% A %*% N0 
A %*% A %*% A %*% N0 
# Det samme beregnet ved hjelp av diagonalisering
D <- diag(eigen(A)$val)
U <- eigen(A)$vec
U %*% D^3 %*% solve(U) %*% N0 # solve beregner inversen til en matrise

# Lister (sammensatt struktur for data betstående av ulike typer)
minliste <- list(a=1,b="per")     # list lager en liste
minliste$b     # $kompenent eller [[i]] refererer til et element i listen ved navn eller nummer
minliste[[1]]
mintest <- t.test(bmi, mu=22.5)     # Funskjoner for ulike statistiske metoder returner typisk en liste
mintest$statistic     
mintest$p.value

# Data frames (en slags matrise og liste nyttig for å representere innsamlede data).
data <- read.table("http://www.math.ntnu.no/~jarlet/statmod/heights.dat")     # Leses typisk inn fra en ekstern fil med read.table eller read.csv
data$sex     # både $komponent og [i,j] kan brukes for å refere til enkeltelement
data[,2]
lm(data$offspring ~ data$mother)
attach(data)     # attach (merk eventuelle advarsler)
lm(offspring ~ mother)
detach(data)     # detach
# search viser hvilke dataframes (og pakker etc) som er attached
 
# Indeksering (hvordan referere til et utvalg element i en vektor, matrise, liste eller dataframe) 
bmi[2]    # ved enkelttall
bmi[c(2,3,2)]    # ved en numerisk vektor for å velge flere element (ett eller flere ganger)
bmi[-c(2,3)]    # negative indekser ekskludere element
vekt[bmi > 22.5]    # ved logiske vektorer laget med logiske operatorer >,<,>=,<=,!=,==
data[ data$sex==1 , ]
subset(data,sex==1)    # (kan også gjøres med subset funksjonen)
minmodell$coefficients["høyde"] # ved navn



# d-, p-, q-, og r- funksjoner
par(mfrow=c(2,2)) # deler opp grafikkvinduet i 2x2 = 4 plot
curve(dnorm(x, mean=2,sd=1.5),-5,10) # tetthetsfunksjonen
curve(pnorm(x, mean=2,sd=1.5),-5,10) # kumulativ fordeling
curve(qnorm(x, mean=2,sd=1.5),0,1) # kvantiler og kvantilfunksjonen
hist(rnorm(n=10000,mean=2,sd=1.5)) # pseudo-tilfeldige tall

par(mfrow=c(2,2)) # deler opp grafikkvinduet i 2x2 = 4 plot
curve(dbinom(x, size=10,prob=0.3),0,10,type="h")
curve(pbinom(x, size=10,prob=0.3),0,10,type="s")
curve(qbinom(x, size=10,prob=0.3),0,1,type="s")
hist(rbinom(n=100, size=10,prob=0.3), 0:11-.5) # pseudo-tilfeldige tall


# set.seed kan brukes for å generere de samme pseudotilfeldige tallene på nytt
set.seed(1)
rnorm(10, mean=2, sd=1.5)

# nytteverdi: finne egenskapene til estimator, gjennomføre numeriske eksperimenter
alpha <- 4 
beta <- 0.5
sigma <- 1
x <- 1:10
betahat <- c()
for (i in 1:1000) {
  Y <- rnorm(10, mean=  alpha + beta*x   ,sd=sigma)
  modell <- lm(Y ~ x)
  betahat[i] <- modell$coef[2]
}
hist(betahat)

# numerisk over- og underflow (problem 3, øving 1)
factorial(365)
factorial(100)
.Machine$double.xmax #(største dobbelpresisjon flyttall)
.Machine$double.xmin #(minste dobbelpresisjon flyttall)

##                                   Mandag 14. jan

# prediksjons- og konfidensbånd
library(ISwR) # Les inn R-pakke som følger læreboka.
thuesen
attach(thuesen)
plot(blood.glucose, short.velocity, ylim=c(.5,2))
# hvordan påvirker glukose-nivå i blodet hjertefunksjon (short.velocity)?
modell <- lm(short.velocity ~ blood.glucose)
abline(modell)
predict(modell)
nyedata <- data.frame(blood.glucose=seq(0,25,by=1))
nyedata
pred <- predict(modell, newdata=nyedata, inter="conf")
lines(nyedata$blood.glucose, pred[,2])
lines(nyedata$blood.glucose, pred[,3])
pred <- predict(modell, newdata=nyedata, inter="pred")
lines(nyedata$blood.glucose, pred[,2])
lines(nyedata$blood.glucose, pred[,3])

detach(thuesen)


## Illustration of the the behaviour of the estimated regression model and 
## confidence interval for a future expected response by simulation
## (You do not need to understand the following code...)
par(mfrow=c(2,1))
n <- 0
hits <- 0
beta <- NULL
repeat
{
  x <- 1:10
  y <- rnorm(n = 10, mean = 5 + 0.25*x, sd = 1)
  plot(x,y,col="red",ylim=c(2,12),xlim=c(-2,13))
  abline(5,.25,col="darkblue") # sann regresjonslinje
  estregr <- lm(y~x)
  abline(estregr,col="red") # estimert regresjonslinje basert på simulerte data
  xnew <- 11; ynew <- 5+.25*xnew
  ci <- predict(estregr,newdata=data.frame(x=xnew),int="c")[2:3]
  lines(rep(xnew,2),ci,col="red")
  n <- n + 1
  beta[n] <- coef(estregr)[2]  
  hist(beta,xlim=c(-.5,1),breaks=seq(-.5,1,by=.05),prob=TRUE,ylim=c(0,7))
  varbeta <- 1/sum((x-mean(x))^2)
  curve(dnorm(x,mean=0.25,sd=sqrt(varbeta)),add=TRUE)
  if (ci[1]ynew)
    hits <- hits + 1
  cat("n = ",n,", hits = ",hits,", hits/n=",hits/n,"\n")
  Sys.sleep(.5) # en liten kunstpause
}

##                          Torsdag 21

# Plotting av funksjoner og parameteriske kurver og punktmengder
t <- seq( 0, 20*pi, length=500)
x <- t*cos(t)
y <- t*sin(t)
plot(x,y,type="l")

## Multiple regression

## plot3d method for lm-objects (fitted multiple regression
## models with two covariates).  Don't try to understand the code
## unless you're feeling adventureous...  Used for pedagogical
## purposes during the lectures
library(rgl)
plot3d.lm <- function(lmobj,plane=TRUE,residualcolor="blue",...) {
  mf <- model.frame(lmobj)
  plot3d(mf[,2],mf[,3],mf[,1],box=F,axes=F,
         xlab=names(mf)[2],ylab=names(mf)[3],zlab=names(mf)[1],...)
  axes3d(c('x','y','z'))
  if (plane) {
    xmin <- min(mf[,2])
    xmax <- max(mf[,2])
    ymin <- min(mf[,3])
    ymax <- max(mf[,3])
    newdata <- data.frame(c(xmin,xmax,xmax,xmin),c(ymin,ymin,ymax,ymax))
    names(newdata) <- names(mf)[2:3]
    ypred <- predict(lmobj,newdata=newdata)
    quads3d(newdata[,1],newdata[,2],ypred,col="white")
    points3d(mf[,2],mf[,3],fitted(lmobj),col=residualcolor)
    for (i in 1:nrow(mf))
      lines3d(rep(mf[i,2],2),rep(mf[i,3],2),c(mf[i,1],fitted(lmobj)[i]),col=residualcolor)
  }
}

##                                Mandag 25. jan.

## Multiple regression, intro, interpretation of the model summary
library(ISwR)
data(cystfibr)
cystfibr
attach(cystfibr)
multregr <- lm(pemax ~ age + bmp)
plot3d(multregr)
summary(multregr)

## Two categories encoded using one dummy variable
multregr <- lm(pemax ~ age + sex)
plot3d(multregr)
summary(multregr)


## Fisher's F-fordeling
par(mfrow=c(2,2))
curve(df(x,df1=10,df2=20),0,5) # tetthetsfunksjonen 
curve(pf(x,df1=10,df2=20),0,5) # kumulativ fordeling
u <- rchisq(10000,df=10)
v <- rchisq(10000,df=20)
hist( (u/10)/(v/20))
hist(rf(10000,10,20)) # det samme...

## To utvalg: Like varianser?
heights <- read.table("http://www.math.ntnu.no/~jarlet/statmod/heights.dat")
heights
attach(heights)
x <- offspring[sex==1]
y <- offspring[sex==0]
var(x)
var(y)
# Testen utført manuelt
f <- var(x)/var(y)
f
qf(c(0.025, 0.975), df1=length(x)-1, df2=length(y)-1) # kritiske verdier
2*pf(f,df1=length(x)-1, df2=length(y)-1, lower.tail=FALSE) # p-value
# var.test
var.test(offspring[sex==1],offspring[sex==0]) # the samme v.h.a. var.test

## One-way analysis of variance
## Example 1
library(ISwR)
data(red.cell.folate)
attach(red.cell.folate)
stripchart(folate~ventilation,vert=T)
## En-veis variansanalyse gjort "manuelt"...
grandmean <- mean(folate) ## total utvalgsgjsnitt
grandmean
groupmean <- tapply(folate,ventilation,mean)    # gruppegjennomsnitt
groupmean
ngroup <- tapply(folate,ventilation,length)  # antall i hver gruppe
ngroup
SSD_T <- sum((folate-grandmean)^2)
SSD_B <- sum(ngroup*(groupmean-grandmean)^2)
SSD_W <- sum((folate - groupmean[ventilation])^2)
SSD_T
SSD_B
SSD_W
MS_T <- SSD_T/(22-1)
MS_B <- SSD_B/(3-1)
MS_W <- SSD_W/(22-3)
MS_T
MS_B
MS_W
Fobs <- MS_B/MS_W
Fobs
qf(0.05, df1=2, df2=19, lower.tail=FALSE)
curve(df(x,df1=2,df=19),0,10)  # a plot of the F-distribution under H_0
abline(v=MS_B/MS_W) # the observed F-value
## Fitting the model and producing the analysis of variance table
folmodl <- lm(folate ~ ventilation)
anova(folmodl)
detach(red.cell.folate)

## En- og to-veis variansanalyse
library(ISwR)
heart.rate
attach(heart.rate)
table(subj,time)
subj
time
mod <- lm(hr ~ time)
anova(mod)
summary(mod)
model.matrix(~time)
interaction.plot(time, subj, hr)
mod <- lm(hr ~ time + subj)
anova(mod)
summary(mod)


## Ikke-balansert design, kolinearitet
x1 <- rnorm(50)
x2 <- .8*x1 + rnorm(50,sd=0.2)
plot(x1,x2)
cor(x1,x2)
y <- 0.5*x1 + 0*x2 + rnorm(50,sd=1)
kolin <- lm(y ~ x1 + x2)
summary(kolin)
plot3d(kolin)
summary(lm(y ~ x1)) # betahat_1 changes when x2 is left out
## vs. balansert design
x1 <- rep(1:5,times=5)
x2 <- rep(1:5,each=5)
table(x1,x2)
cor(x1,x2) 
y <- 0.5*x1 + 0*x2 + rnorm(25,sd=1) # same model 
balans <- lm(y~x1+x2)
plot3d(balans)
summary(balans)
summary(lm(y~x1)) # betahat_1 remain unchanged when x2 is left out

##                          Torsdag 4. februar

## Model seleksjon basert på F-tester
set.seed(3)
f1 <- factor(sample(c("a","b","c","d"),20,repl=T)) # a factor
x1 <- rnorm(20,(1:4)[f1])
x2 <- rnorm(20,x1,sd=.5)
y <- rnorm(20,
           mean=10 + 0.01*x1 + 0.2*x2 + c(0,1,1.5,-.5)[f1],
           sd=.5)
completedata <- data.frame(y,x1,x2,f1)
rm(x1,x2,f1)
trainingset <- completedata[1:10,]
validationset <- completedata[11:20,]

trainingset
fullmodel <- lm(y ~ x1 + x2 + f1,data=trainingset)
summary(fullmodel) 
# forenklet modell
reduced1 <- lm(y ~ x1 + x2, data=trainingset)
# F-test v.h.a. anova
anova(reduced1, fullmodel)
# drop1 gir F-tester av alle innskrekninger mot aktuell modell
drop1(fullmodel, test="F")   
reduced1 <- lm(y ~ x2 + f1, data=trainingset)
drop1(reduced1, test="F")
reduced2 <- lm(y ~ f1, data=trainingset)   
drop1(reduced2, test="F")
# add1 gir F-tester av aktuell model mot ulike utvidelser angitt i scope
add1(reduced2, .~.+x1+x2, test="F")
drop1(lm(y ~ x1, data=trainingset),test="F")
drop1(lm(y ~ x2, data=trainingset),test="F")

step(fullmodel)


# Hvor gode prediksjoner gir valgt modell sammenlignet med den fulle modellen?
# Er enkle modeller bedre enn kompliserte modeller?
predict(fullmodel, newdata=validationset)   # prediksjoner basert på full modell
predict(reduced, newdata=validationset)   # prediksjoner basert på valgt modell
validationset$y   # de reelle observasjonene
# Mean square prediction error for fullmodel   
mean((validationset$y - predict(fullmodel, newdata=validationset))^2)
# Mean square prediction error for reduced2 (valgt model)
mean((validationset$y - predict(reduced2, newdata=validationset))^2)

# F-testene vi får om vi kjører anova på modell uten balansert design avhenger av
# rekkefølge og er skjeldent av interesse
anova(fullmodel)
anova(lm(y ~ f1 + x1 + x2, data=trainingset))

## AIC



## Flere modeller kan være virke rimelige (tilbake til kolinearitetseksempel).  
x1 <- rnorm(50)
x2 <- .8*x1 + rnorm(50,sd=0.2)
plot(x1,x2)
cor(x1,x2)
y <- 0.5*x1 + 0*x2 + rnorm(50,sd=1)
fullmodel <- lm(y ~ x1 + x2)
drop1(fullmodel, test="F")
summary(fullmodel)

drop1(lm(y ~ x1), test="F")
drop1(lm(y ~ x2), test="F")


# Multinomisk fordeling:
#   Plot av simultanfordelingen P(X1=x1 and X2=x2) og 
#   marginalfordelingene til X1 og X2 når
#   X1,X2,X3 ~ multinom(.2, .6, .2, 10)
#   (koden under er uviktig)
n <- 20 # endre til 100
p <- c(.1,.6,.3)
layout(mat = matrix(c(2,1,0,3),2,2,byrow=TRUE),widths = c(3,5), heights = c(5,3))
plot(NA,xlim=c(0,n),ylim=c(0,n),xlab=expression(x[1]),ylab=expression(x[2]))
for (x1 in 0:n) 
  for (x2 in 0:(n-x1)) 
    points(x1,x2,cex=15*sqrt(dmultinom(c(x1,x2,n-x1-x2),size=n,prob=p)))
text(6,6,expression(P(X[1]==x[1],X[2]==x[2])))
plot(NA,xlim=c(0,max(px2)),ylim=c(0,n),xlab=expression(P(X[2]==x[2])),ylab=expression(x[2]))
segments(0,0:n,dbinom(0:n,size=n,prob=p[2]),0:n)
plot(0:n,dbinom(0:n,size=n,prob=p[1]),type="h",xlab=expression(x[1]),ylab=expression(P(X[1]==x[1])))

## Test av Benfords lov 
data <- read.table("http://www.math.ntnu.no/~jarlet/statmod/mammals.dat",header=TRUE)
X <- data$body
X
X <- floor(X/10^floor(log10(X))) # divider på nærmeste heltallige 10-potens og rund av nedover
X
X <- table(factor(X,levels=1:9)) # tell opp antall 1'ere, 2'ere, ... 9'ere.
X # observerte antall
n <- sum(X) # totalt antall delforsøk i den multinomiske forsøksrekka
i <- 1:9
p <- log(i+1, 10) - log(i, 10) # teoretiske sannsynligheter
p
n*p # forventet antall av hvert siffer
D <- sum((X - n*p)^2/(n*p)) # testobservator
D
qchisq(.05,df=8,lower.tail=FALSE) # kritisk verdi
pchisq(D,df=8,lower.tail=FALSE)
chisq.test(X,p=p) # det samme med chisq.test()
# uniform fordeling kan forkastes
chisq.test(X,p=rep(1/9,9))

## Test av uavhengighet i 2 x 2 kontingenstabell
X <- matrix(c(23,18,7,13),2,2)
X
chisq.test(X, correct=FALSE)
chisq.test(X, simulate.p.value = TRUE,B=1e+5) # Mere eksakt
chisq.test(X, correct=FALSE)$observed
chisq.test(X, correct=FALSE)$expected

## Test av uavhengighet i r x c kontingenstabell
library(ISwR)
?juul
attach(juul)
table(tanner,sex)
chisq.test(tanner,sex)$observed
chisq.test(tanner,sex)$expected

# Testing H-W equilibrium
X <- c(51,42,7) # number of AA, Aa, aa in the sample
n <- sum(X)
phat <- (2*X[1] + X[2])/(2*n)
phat
Phat <- c(phat^2, 2*phat*(1-phat), (1-phat)^2)
Phat
n*Phat
D <- sum((X - n*Phat)^2/(n*Phat))
D

#             Mandag 15. feb.

# Incomplete data, blood type example 
x <- c(44,27,4,88) # Fra Crow, Introduction to population genetics theory
n <- sum(x)
multinomialprobs <- function(par) {
  pA <- par[1]
  pB <- par[2]
  p0 <- 1-pA-pB
  c(pA^2 + 2*pA*p0,  # freq. of bloodtype A
    pB^2 + 2*pB*p0,  # freq. of bloodtype B
    2*pA*pB,         # freq. of bloodtype AB
    p0^2)            # freq. of blootype 0
}
n*multinomialprobs(c(.2,.3))
lnL <- function(par,x) {
  n <- sum(x)
  -dmultinom(x,prob=multinomialprobs(par),log=T)
}
fit <- optim(c(.25,.25),lnL,x=x)
fit
fit$par # estimated frequency of allele A and B
phat <- multinomialprobs(fit$par)
phat # estimated phenotype frequencies
D <- sum((x-phat*n)^2/(n*phat))
D # goodness of fit statistics
pchisq(D,df=4-1-2,lower.tail=FALSE) # p-value


# Fisher's eksakte vs andre alternativ
X <- matrix(c(20,10,10,0),2,2)
X
chisq.test(X) # Multinomisk fordeling og kji-kvadrat tilnærming
chisq.test(X,simulate.p.value=TRUE,B=1e+4) # Multinomisk og simulering 
fisher.test(X) # Fishers test basert på konstante marginaler


# Generaliserte lineære modeller
# simulated example of logistic regression with one covariate
x <- 0:10
n <- length(x)
eta <- -2.5 + 0.8*x
p <- 1/(1+exp(-eta))
y <- rbinom(n=n, prob=p, size=10)
plot(x,y/10,ylim=c(0,1))
curve(1/(1 + exp(-(-2.5 + 0.5*x))), add=TRUE)


#                 Torsdag 17. febr.
        
# Dalg. p. 229
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)
data.frame(smoking,obesity,snoring,n.tot,n.hyp)
# n x 2 matrix containing number of success and failures on the left hand side
n.tab <- cbind(n.hyp, n.tot-n.hyp)
n.tab
logreg3 <- glm(n.tab ~ smoking + obesity + snoring, fam=binomial())
summary(logreg3)
logreg2 <- glm(n.tab ~ obesity + snoring, fam=binomial())
summary(logreg2)
deviance(logreg2) - deviance(logreg3)
# kritisk verdi
qchisq(0.05, df=4-3, lower.tail=FALSE)
# drop1 med kji-kvadrattest
drop1(logreg3, test="Chisq")


#                 Mandag 22. febr.
 
## Mettet (saturated) modell
set.seed(1)
n <- 10 # antall delforsøk i hver observasjon
y <- rbinom(n=5,size=n,prob=.5)
x <- 1:5
plot(x,y/n,ylim=c(0,1),xlim=c(0,6))
yy <- cbind(y, 10-y)
xseq <- seq(0,6,length=100)
saturated <- glm(yy ~ x + I(x^2) + I(x^3) + I(x^4), binomial)
lines(xseq,predict(saturated, newdata=data.frame(x=xseq),type="response"),col="red")
predict(saturated, type="response")
y/n

## Dalg. p. 239
library(ISwR)
data(juul)
juul$menarche <- factor(juul$menarche, labels=c("No","Yes"))
juul$menarche
juul.girl <- subset(juul, age>8 & age< 20 & complete.cases(menarche))
attach(juul.girl)
plot(age,menarche=="Yes")
menmod <- glm(menarche ~ age,binomial)
summary(menmod)
beta <- coef() # SME av intercept og stigningstall beta1 og beta2
curve(1/(1+exp(-(beta[1]+beta[2]*x))),add=TRUE)

## standardfeil til x0hat (gjennomsnittlig alder ved menarche)
vc <- vcov(menmod) # varians-kovariansmatrise
vc
der <- c(-1/beta[2],beta[1]/beta[2]^2 ) # partiellderiverte til x0 = f(beta1, beta2)
# varians til x0hat
varx0hat <- (der[1])^2*vc[1,1] + (der[2])^2*vc[2,2] + 2*der[1]*der[2]*vc[1,2]
t(der) %*% vc %*% der # det samme uttrykkt i matrisenotasjon
# og standard feil
sqrt(varxh0hat)

#                  Torsdag 25. febr

# probit link (situasjon 1)
menmod <- glm(menarche ~ age,binomial(link="probit"))
summary(menmod)
# Omregnet til mu og sigma får vi:
1/0.86 # sigma
-(-11.37)/0.86


# probit link (situasjon 2)
# simulert eksempel
par(mfrow=c(2,1))
x <- runif(1000)
beta0 <- 1
beta1 <- 4
sigma <- 1.5
Z <- rnorm(1000, mean = beta0 + beta1*x, sd=sigma)
plot(x,Z)


threshold <- 2
abline(h=threshold)
Y <- ifelse(Z > threshold, 1, 0)
plot(x,Y+runif(1000,0,.1),pch=".")
summary(glm(Y ~ x, binomial(link="probit")))
(beta0 - threshold)/sigma
beta1/sigma


#                 Mandag 29. febr

## cloglog eksempel
sparrows <- read.table("http://www.math.ntnu.no/~jarlet/statmod/cloglog_data_sparrows_Plos_ONE.txt",header=TRUE)
sparrows
attach(sparrows)
propdead <- dead/(dead+alive)
plot(para_eggs,propdead)
plot(hatchday,propdead)
plot(chickage,propdead)
yy <- cbind(dead,alive)
summary(glm(yy ~ para_eggs + hatchday + log(chickage), binomial(link="cloglog")))
# Estimat av dødsrate (dimensjon dag^-1)
exp(0.073 + 0.007*0 + -0.0186*174)
detach(sparrows)


## Poisson regression, eksempel 1
library(ISwR)
eba1977
attach(eba1977)
fit <- glm(cases ~ city + age, offset=log(pop), fam=poisson)
summary(fit)
## Interpretation of parameter estimates

## Goodness-of-fit reasonable (cannot reject the model)

## Which factors should be included? (drop1)
drop1(fit, test="Chisq")
## Fit the model assuming that only Fredericia is different due the petrochemical factory (city=="Fredericia")
frederica <- factor(city=="Fredericia")
frederica
fit2 <- glm(cases ~ frederica + age, offset=log(pop), fam=poisson)
summary(fit2)
drop1(fit2,test="Chisq")
anova(fit,fit2,test="Chisq")
# We conclude that fit2 is the best model


#                    Torsdag 3. mars

## Poisson regresjon, eksamen 2015v, oppgave 2, Hoge Veluwe data fra 2008-sesongen
load(url("https://www.math.ntnu.no/~jarlet/statmod/hv2008.RData"))
attach(hv2008)
plot(t,y)
t2 <- t^2
mod <- glm(y ~ t + t2, family=poisson(link="log"))
summary(mod)
beta <- coef(mod)
curve(exp(beta[1] + beta[2]*t + beta[3]*t^2), xname="t",add=TRUE)
qchisq(0.05, df=89, lower.tail=FALSE) # deviance larger than critical value suggests over-dispersion
pchisq(deviance(mod),df=df.residual(mod),lower.tail=FALSE)

#                    Mandag 7. mars

# Korreksjon for overdispersjon
# Estimat av dispersjonsparameter beregnet "manuelt":
sum(residuals.glm(mod,type="pearson")^2)/df.residual(mod)
# Det samme beregnet ved å bruke quasipoisson
mod2 <- glm(y ~ t + t2, family=quasipoisson(link="log"))
summary(mod2)
# Sammenligning av modeller, F-tester i stedet for kji-kvadrat
# tester basert på endring i devians:
drop1(mod, test="F")
detach(hv2008)



## Lineær separasjon
## Eksempel 1 (fra øving 7)
data <- read.csv("~/Downloads/wildlife-encounters - Sheet1 (3).csv")
attach(data)
summary(glm(moose ~  sex + studyprogram + hours, offset=log(t), family=binomial(link="cloglog")))
detach(data)
## Eksempel 2
x <- 1:10
y <- rep(c(0,1),each=5)
plot(x,y)
summary(glm(y ~ x,binomial))
## Eksempel 3, se notat 4.


## interactions between one factor and one numeric covariate
library(ISwR)
hellung
attach(hellung)
names(hellung)
plot(conc,diameter,col=glucose)
plot(conc,diameter,col=glucose+1,log="xy")
reg1 <- lm(log10(diameter) ~ log10(conc), data=hellung[glucose==1,])
reg2 <- lm(log10(diameter) ~ log10(conc), data=hellung[glucose==2,])
abline(reg1)
abline(reg2)
summary(reg1)
summary(reg2)
## approximate test based on normality for difference in slope comparing estimates from each model fit
z <- (-.0596 + 0.05320)/sqrt(0.004125^2 + 0.00272^2)
z
# med interaksjon
fglucose <- factor(glucose)
fglucose
samspill <- lm(log10(diameter) ~ log10(conc) + fglucose + log10(conc):fglucose, data=hellung)
summary(samspill)
drop1(samspill,test="F")
detach(hellung)


## interaksjon mellom to (eller flere) faktorer
heli <- read.csv("~/Downloads/helicopterdata - Sheet1.csv")
attach(heli)
mod.y <- lm(flighttime ~ wing + clip + wing:clip)
anova(mod.y)
summary(mod.y)
par(mfrow=c(2,1))
interaction.plot(wing,clip,flighttime,ylim=c(0,15),lty=1,col=1:2)
points(wing,flighttime,col=clip)

mod.lny <- lm(log(flighttime) ~ wing + clip)
anova(mod.lny)
summary(mod.lny)
interaction.plot(wing,clip,log(flighttime),lty=1,col=1:2,ylim=c(1.3,2.5))
points(wing,log(flighttime),col=clip)



## some simple programming examples
## -------------------------------------------------------------------------
## weighted average
x <- c(1,2,3,4,5)
mean(x) # ordinary arithmetic mean
w <- c(.4,.3,.2,.1,0)  # weights
sum(w)  
sum(w*x)
# weighted mean

## wmean:
## Description: Computed a weighted average
## 
## Arguments:
##   x: A numeric vector
##   w: A numeric vector containing the weights
## Value:
##   The weighted mean of x
## Examples:
wmean <- function(x, w) {
  sum(w*x)
}
## Tests:
wmean(c(1,2,3,4,5), c(.4,.3,.2,.1,0))


## -------------------------------------------------------------------------
## quad.eq.solve 
## Description: Computes the solution of a quadratic equation
##
##   a x^2 + b x + c = 0
##
## (illustrate use of if-statements and control-flow)
##
## Arguments:
##   a,b,c : The coefficient a, b, c
## Value:
##   A vector with the real solutions as elements
## Examples:
##   quad.eq.solve(1,-2,1) -> 1  
##   quad.eq.solve(1,-2,3) -> Empty vector
##   quad.eq.solve(1,-3,2) -> A vector with elements 1 and 2
quad.eq.solve <- function(a,b,c) {
  delta <- b^2 - 4*a*c
  if (delta>0) {
    sol <- c((-b + sqrt(delta))/(2*a),  (-b - sqrt(delta))/(2*a) )
  } else { # delta <= 0
    if (delta==0) {  
       sol <- -b/(2*a)
    } else { # delta < 0, no solutions
      sol <- c()  
    }
  }
  sol  
}
## Test that it works (by substitution)
quad.eq.solve(1,-2,1)
quad.eq.solve(1,-2,3)
x <- quad.eq.solve(1,-3,2)
x^2 - 3*x + 2

## For-loops.  Compute the first 10 numbers of the Fibbonacci sequence
x <- c(1,1)
for (i in 3:10) {
  x[i] <- x[i-1] + x[i-2] 
}
x

## -------------------------------------------------------------------------
## fib
## Description: A function returning the first n numbers of the Fibbonacci sequence
##
## Arguments:
##
## Value:
##
## Examples:
##


## Tests:

## -------------------------------------------------------------------------
## simulate a deterministic discrete time logistic population dynamic model
N <- 10
K <- 100
r <- 0.1
for (t in 2:100)
  N[t] <- N[t-1] + r*N[t-1]*(1-N[t-1]/K)
plot(N,type="l")

## logistic growth
## Description: Simulate a discrete logistic population dynamic model
##
## Arguments:
##   N1: Initiell populasjonstørrelse
##   r: vektstrate
##   K: bæreevne
##   tmax: Antall generasjoner
## Value:
##   En vektor som inneholder N_1, N_2, ... N_tmax
logisticgrowth <- function(N1, r, K, tmax) {
  N <- rep(NA,tmax)
  N[1] <- N1
  for (t in 2:tmax) {
    N[t] <- N[t-1] + r*N[t-1]*(1-N[t-1]/K)
  }
  N
}
## Tests:
plot(logisticgrowth(10, 0.2, 1000, 200))



##                  Torsdag 31. mars


## -------------------------------------------------------------------------
## Mysqrt
## Description: Find the square root of a using Newton's method (illustrate
## while-loops)
## 
## Arguments:
##   a: 
## Value:
##   The square root of a
## Examples:
##   mysqrt(2) -> 1.414
##   mysqrt(4) -> 2
mysqrt <- function(a) {
   x <- a/2
   while (abs(x^2 - a) > 1e-8) {
      x <- x/2 + a/(2*x)
      print(x)
   }  
   x  
}
## Tests:
mysqrt(2)




##                         Mandag 4. april

## Programming and numerical maximisation of the likelihood function
tsim <- rweibull(200,shape=2,scale=10)
hist(tsim,breaks=10,freq=F)
curve(dweibull(x,shape=2,scale=10),add=T)

## log likelihood funksjon for Weibull-fordelte data
##
## Argument
##   par: En vektor som inneholder verdien av skala og form-parameter
##   t: En vektor som innneholder dataene t_1,..., t_n
## Verdi
##   Verdien av (det negative) log likelihoodet gitt a, b og t'ene.
lnL <- function(par,t) {
  a <- par[1]
  b <- par[2]
  n <- length(t)
  -(n*log(a) -n*log(b) + (a-1)*sum(log(t/b)) - sum((t/b)^a))
}
## Enklere versjon som bruker dweibull
lnL <- function(par,t) {
  a <- par[1]
  b <- par[2]
  -sum(dweibull(t, shape=a, scale=b, log=TRUE))
}

## Numeriks maksimering av log likelihood funksjonen
optim(c(1,.1), lnL, t=tsim)


## Overflateplot av likelihood-funksjonen (motivasjon av standardfeil basert på asymptotisk teori)
library(emdbook)
## Hva er effekten av økende utvalgsstørrelse?
tsim <- rweibull(25,shape=2,scale=10)
curve3d(exp(-lnL(c(x,y),t=tsim)),
        xlim=c(1,3),ylim=c(5,15),        
        sys="persp",
        theta=30,phi=30,
        tick="detailed",
        xlab="shape parameter a",ylab="scale parameter b",zlab="likelihood")


# Tilnærmede standard feil based on asymptotic theory
fit <- optim(c(1,.1), lnL, t=tsim, hessian=TRUE)
fit
sqrt(diag(solve(fit$hessian)))
fit$par

# Test av nøstede modeller basert på 2 ganger endring i maksimalt log likelihood
lnL1 <- -fit$value
n <- length(tsim)
lnL0 <- -n*(log(mean(tsim))+1)
lnL0
lnL1
2*(lnL1-lnL0)
qchisq(.05, df=1, lower.tail=FALSE)



##                     Torsdag 7. april

## Multinomiske blodtypedata, goodness of fit
x <- c(44,27,4,88)
n <- sum(x)
n

## Likelihoodfunksjon for multinomiske data
##
## Argument:
##    par: En vektor som inneholder allelfrekvensene p_A, p_B eller tilsvarende
##    x: En vektor som inneholder dataene x_1,...,x_k
##    probfn: En funksjon som beregner de multinomiske sannsynlightene fra et mindre
##       antall parametere inneholdt i vektoren par 
## Verdi:
##    (det negative) log likelihoodet
lnL <- function(par, x, probfn) {
  p <- probfn(par)
  -dmultinom(x,prob=p,log=TRUE)
  # Alternativt:
  # n <- sum(x)
  # -(lfactorial(n) - sum(lfactorial(x)) + sum(x*log(p)))
}
## Fenotypefrekvensene som funksjon av allelfrekvensene for hypotese 1
## Argument:
##   par : En vektor som inneholder allelfrekvensen p_A og p_B
## Verdi:
##   En vektor som inneholder de 4 fenotypefrekvensen
multinomialprobs <- function(par) {
  pA <- par[1]
  pB <- par[2]
  p0 <- 1-pA-pB
  c(pA^2 + 2*pA*p0, pB^2 + 2*pB*p0 , 2*pA*pB , p0^2)
}

fit <- optim(c(.1,.1),lnL, x=x, probfn=multinomialprobs, hessian=TRUE)
fit$par # SME av p_A og p_B
sqrt(diag(solve(fit$hessian))) # Standardfeil til p_A og p_B
cov2cor(solve(fit$hessian)) # Korrelasjon mellom SME av p_A og p_B
# observed value of chi square statistic
phat <- multinomialprobs(fit$par) # SME av frekvens av observarbare fenotyper (blodtypene)
n <- sum(x)
D <- sum((n*phat - x)^2/(n*phat))
D

## Generalisert lineær modell eksempel
library(ISwR)
data(juul)
juul$menarche <- juul$menarche - 1
juul.girl <- subset(juul, age>8 & age< 20 & complete.cases(menarche))
attach(juul.girl)
plot(age,menarche)
# SMEer og standardfeil beregnet av glm
summary(glm(menarche ~ age,fam=binomial)) 
curve(1/(1+exp(-(-20.01+1.5173*x))),add=T)

## Samme modell tilpasset v.h.a. hjemmelaget likelihoodfunksjon og
## numerisk maksimalisering av denne.  Merk at dette gir eksakt samme resulatet
## både når det gjelder estimater og deres standardfeil

# Likelihood funksjon for generalisert lineær modell med
# én numerisk forklaringsvariabel (utgangspunkt for problem 2, øving 10)
#
# Argument
#   par: En vektor som inneholder verdiene på beta_0 og beta_1
#   y: En vektor som inneholder response y_1,y_2,...,y_n
#   x: En vektor som inneholder verdien på forklaringsvariabel x_1,...,x_n
# Verdi
#   (det negative) log likelihoodet
lnLglm <- function(par, y, x) {
  beta0 <- par[1]
  beta1 <- par[2]
  p <- 1/(1 + exp(-(beta0 + beta1*x)))
  -sum(dbinom(y, size=1, prob=p, log=TRUE))
  # Alternativt -sum(y*log(p) + (1-y)*log(1-p)) men dette kan numerisk gi log(0)=-Inf
  # hvis noen av p_i'ene er (numerisk) lik 0
}
fit <- optim(c(0,0), lnLglm, hessian=TRUE, x=age, y=menarche)
fit$par # SME av intercept og slope
sqrt(diag(solve(fit$hessian))) # tilhørende standardfeil


#                       Mandag 11. april

# Maksimalisering av likelihoodet med begrensninger på parameterne
optim(c(.5,.1),lnL,x=x,probfn=multinomialprobs)

# Ved reparameterisering: Modifisert programkode hypotese 1, øving 10
# Lar par[1] og par[2] svare til theta[1] og theta[2]
multinomialprobs <- function(par) {
  pA <- exp(par[1])/(1+exp(par[1])+exp(par[2])) 
  pB <- exp(par[2])/(1+exp(par[1])+exp(par[2]))
  p0 <- 1-pA-pB
  c(pA^2 + 2*pA*p0, pB^2 + 2*pB*p0 , 2*pA*pB , p0^2)
}
optim(c(10,1),lnL,x=x,probfn=multinomialprobs)





 
#
#     Beregning the maksimalt log likelihood under mettet model
#
# Eksempel, n=4 observasjon, p = 1, 2, 3, 4 parametere.
n <- c(10,20,15,1)
x <- -1:2
set.seed(2)
y <- rbinom(4,size=n,prob=plogis(.1+.1*x1))
plot(x,y/n,cex=sqrt(n),ylim=c(0,1))
y/n
## fitted values for models with 1, 2, 3 and 4 parameters:
yy <- cbind(y,n-y)
lines(x1,predict(glm(yy ~ 1, binomial),type="response"))
lines(x1,predict(glm(yy ~ x, binomial),type="response"))
lines(x1,predict(glm(yy ~ x + I(x^2), binomial),type="response"),col="red")
lines(x1,predict(metted <- glm(yy ~ x + I(x^2) + I(x^3), binomial),type="response"),col="blue")
metted
logLik(metted)

# Reparameterisering med p_1, p_2, ..., p_3 gir SMEer lik y_i/n_I og lnL_saturated gitt
# ved (det ikke helt generelle uttrykket)
sum(lchoose(n,y) + y*log(y/n) + (n-y)*log(1-y/n))
# eller, for å håndere at log(0^0)=1 og ikke NaN,
sum(lchoose(n,y) + log((y/n)^y) + log((1-y/n)^(n-y)))
# eventuellt
sum(dbinom(y, size=n, prob=y/n, log=TRUE))



#                  Torsdag 14. april

# Forventningsverdier og sannsynligheter beregnet ved simulering
# som grenseverdi til n_A/n og 1/n sum(X) når n -> uendelig
cummean <- function(x) cumsum(x)/(1:length(x))
plot(cummean(rnorm(1000,mean=5,sd=2)),type="l")    # rnorm
# Hvis U ~ unif(0,1), hva er da E(1/U)?
u <- runif(100000)
plot(cummean(1/u),type="l")

# (1-alpha)-Konfidensintervall (A,B) for theta. Er P(A < theta < B) = 1 - alpha?
# Anta at X ~ bin(n,p).  R beregner følgende intervall for p, gitt X=10 og n=100
prop.test(10,100)$conf.int 

# Tilsvarende intervall i Løvås (baserer seg på flere tilnærminger) av pedagogiske hensyn)
løvås <- function(x,n,alpha=0.05) {
  phat <- x/n
  z <- qnorm(alpha/2,lower.tail=FALSE)
  ci <- c(phat - z*sqrt(phat/(1-phat)/n),  phat + z*sqrt(phat/(1-phat)/n))
  list(conf.int=ci)
}
løvås(10,100)

# Dekningsgrad til konfidensintervall
# 
# Argument
#   n: utvalgstørrelse  
#   p: Paramteren p i binomisk modell
#   fn: En funksjon som beregner konfidensintervall i binomisk modell
#   nsim: Antall simulering
# Verdi
#   Estimert deknkingsgrad til konfidensintervall
dekningsgrad <- function(n,p,fn,nsim=1000) {
  nA <- 0
  for (i in 1:nsim) {
    x <- rbinom(1, size=n, prob=p)
    konfint <- fn(x,n)$conf.int
    if (konfint[1] < p & konfint[2] > p) {
       nA <- nA + 1     
    }
  }
  nA/nsim
}
dekningsgrad(100,0.1,fn=prop.test, 1e+4)
dekningsgrad(100,0.1,fn=løvås, 1e+5)


# Varians-kovariansmatrise til SME av a og b for weibullfordelte data estimert 
# ved å simulere (såkalt parameterisk bootstrapping)
lnL <- function(par, t) {
  -sum(dweibull(t,scale=par[1],shape=par[2],log=TRUE))
}
# anta følgende observerte data
t <- rweibull(10,shape=2,scale=10)
hist(t)
# SME av a og b basert på estimerte data
fit <- optim(c(1,1),lnL,lower=c(0,0),t=t,method="L-BFGS-B",hessian=TRUE)
fit$par
# tilnærmet varians-kovarians matrise basert på asymptotisk teori
solve(fit$hessian) 
sqrt(diag(solve(fit$hessian)))
fit$par

# det samme basert på 1000 bootstrapsample
bootstrapestimat <- matrix(NA, 1000, 2)
for (i in 1:1000) {
  # simuler nye "bootstrap"-sample ved pametriske simulering (dweibull)
  bootstrapsample <- rweibull(10, scale=fit$par[1], shape=fit)$par[2])
  # tilpass modellen på nytt og ta vare på bootstrap-estimatene
  bootstrapestimat[i,] <- optim(c(1,1),lnL,lower=c(0,0),t=bootstrapsample,method="L-BFGS-B")$par
}
plot(bootstrapestimat)
points(fit$par[1], fit$par[2], col="red")
sd(bootstrapestimat[,1])
sd(bootstrapestimat[,2])

# ikke-parameterisk bootstrapping
bootstrapestimat <- matrix(NA, 1000, 2)
for (i in 1:1000) {
  # simuler nye "bootstrap"-sample ved å trekke med tilbakelegging fra de orginale dataene
  bootstrapsample <- sample(t,size=50,replace=TRUE)
  # tilpass modellen på nytt og ta vare på bootstrap-estimatene
  bootstrapestimat[i,] <- optim(c(1,1),lnL,lower=c(0,0),t=bootstrapsample,method="L-BFGS-B")$par
}
plot(bootstrapestimat)
points(fit$par[1], fit$par[2], col="red")
sd(bootstrapestimat[,1]) # Standarfeil basert på bootstrapestimatene 
sd(bootstrapestimat[,2])



# Er 2(lnL1-lnL0) kji-kvadrat med p_1 - p_0 frihetsgrader?


#                  Mandag 18. april

# Teststyrke:

# Ett-utvalgs t-test, ikke-sentral t-fordeling
curve(dnorm(x),-5,5)    # standard normalfordeling
curve(dt(x,df=5),-5,5,add=TRUE,col="red")    # t-fording
curve(dt(x,df=5,ncp=2),add=TRUE,col="darkgreen")    # ikke-sentral t-fordeling
curve(pt(x,df=5,ncp=2),-5,5,col="darkgreen")    # kumulativ ikke-sentral t-fordeling
# power.t.test
power.t.test(n=100, delta=.2, sd=1, sig.level=0.05, type="one.sample")
power.t.test(power=0.8, delta=.2, sd=1, sig.level=0.05, type="one.sample")
power.t.test(power=0.8, n=50, sd=1, sig.level=0.05, type="one.sample")

# Proposjoner: Data fra øving 5
tabell <- matrix(c(17,4,11,8),2,2)
colnames(tabell) <- c("left","right")
rownames(tabell) <- c("biology","biotech")
tabell
chisq.test(tabell)
# hvor mye mer data trengs for å påvise en signifikant forskjell (med sannsynlighet 0.8) 
# gitt dataene i "pilot"-studiet?
# power.prop.test
p1 <- 17/(17+11)
p2 <- 4/(4+8)
power.prop.test(p1=p1,p2=p2, power=0.8)

# Hvor robust (vurdert ut i fra P(type I feil)) er en to-utvalgs test hvis dataene i realiteten er
# f.eks. ekponentiellt fordelt med felles parameter lambda under H_0.  Hvordan
# avhenger dette av utvalgsstørrelsen
#
# Argument:
#   n: Felles størrelse av hvert utvalg
#   lambda: Felles parameter i hver ekponentiell fordeling
#   nsim: Antall simulerte utvalg
# Verdi:
#   sannsynlighet for type I feil estimert ved å simulere
robusthet <- function(n, lambda, nsim) {
  n.forkast <- 0
  for (i in 1:nsim) {
     x <- rexp(n,lambda)
     y <- rexp(n,lambda)
     if (abs(t.test(x,y, var.equal=TRUE)$statistic) > qt(.975, df=2*n-2)) {
       n.forkast <- n.forkast + 1
     }
  }
  n.forkast/nsim # estimert P(type I feil)
}
robusthet(100,10,nsim=1e+4)
robusthet(30,10,nsim=1e+4)
robusthet(10,10,nsim=1e+5)
robusthet(3,10,nsim=1e+5)