# KOMMANDOER TIL FORELESNINGENE FOR UKE 4 

# ===============

# EKSEMPEL 11.7 (teststyrke)

# Vi vil se hvordan vi kan gjoere beregningene i eksempel 11.7
# ved hjelp av kommandoen "power.anova.test" (det er ikke 
# noedvendig aa forstaa hvordan en bruker diagrammene paa 
# sidene 574-575 i D&B)

# Setter verdi for antall grupper og antall observasjoer 
# pr gruppe:
I=4
J=8

# Setter verdier for my-ene og sigma:
my=c(0,0,0,1)
s=1

# Bestemmer teststyrken:
power.anova.test(groups=I, n=J, between.var=var(my), within.var=s^2, power=NULL)

# Bestemmer antall observasjoner pr gruppe som gir 95% styrke:
power.anova.test(groups=I, n=NULL, between.var=var(my), within.var=s^2, power=0.95)

# ===============


# EKSEMPEL: OPPTAK AV FETTSYRER
# (fortsettelse fra uke 3)

# Leser inn dataene 
fett=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v16/fettsyrer.txt",header=T) 

# Plott av dataene:
plot(fett$produkt,fett$opptak, ylim=c(0,300),xlab="Produkt",ylab="Opptak")

# Enveis variansanalyse (jf. tabell 11.2 paa side 559 i D&B):
aov.fit=aov(opptak~factor(produkt),data=fett)
summary(aov.fit)

# Simultane konfidensintervall [jf. (11.4) paa side 566 i D&B]:
tukey.fit=TukeyHSD(aov.fit,ordered=T)
print(tukey.fit)
plot(tukey.fit)

# 95% konfidensinterval for kontrsten (my2+my3)/2-(my4+my5)/2
# (jf. nederst side 570 i D&B)
I=5
J=6
snitt=rep(0,I)
for (i in 1:I) snitt[i]=mean(fett$opptak[fett$produkt==i])
cc=c(0,1/2,1/2,-1/2,-1/2)
MSE=sum(aov.fit$residuals^2)/(I*(J-1))
tettahatt=sum(cc*snitt)
se.tettahatt=sqrt(MSE*sum(cc^2/J))
nedre.grense=tettahatt-qt(0.975,I*(J-1))*se.tettahatt
ovre.grense=tettahatt+qt(0.975,I*(J-1))*se.tettahatt
c(nedre.grense,ovre.grense)


# ===============

# EKSEMPEL 11.8 (enveis anova med ulikt antall observasjoner)
# Elstisitetsmodul for tre stoepeprosesser for av magnesium
# https://no.wikipedia.org/wiki/Elastisitetsmodul

# Les inn dataene:
exmp11.8=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v16/exmp11-08.txt", header=T) 

# Boksplott av dataene:
boxplot(ElastMod~Prosess,data=exmp11.8)

# Enveis variananalyse:
fit.exmp11.8=aov(ElastMod~factor(Prosess),data=exmp11.8)
summary(fit.exmp11.8)

# Normalfordelingsplott
qqnorm(fit.exmp11.8$residuals)

# Test for lik varians:
library(car)
leveneTest(ElastMod~factor(Prosess),data=exmp11.8,center=mean)

# Simultane konfidensintervall
tukey.exmp11.8=TukeyHSD(fit.exmp11.8,ordered=T)
print(tukey.exmp11.8)
plot(tukey.exmp11.8)