### Eksempler Kap. 10

# Leser data
rotter = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/rotter.vekt.txt",header=TRUE)
dim(rotter)
names(rotter)
head(rotter)
tail(rotter)
summary(rotter)
summary(rotter$vektok[rotter$beh==1])
summary(rotter$vektok[rotter$beh==2])

# Boxplot
par(mfrow=c(1,1))
boxplot(vektok~beh,data=rotter)

# Histogram
par(mfrow=c(1,2))
hist(rotter$vektok[rotter$beh==1],xlab="vektøkning",nclass=20)
hist(rotter$vektok[rotter$beh==2],xlab="vektøkning",nclass=20)



## Eks. 10.4 - i forrige versjon av Læreboka
# Karakterer i høst og vårsemester

# Leser dataene
score = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp10-05.txt",header=TRUE)

# Boxplot 
boxplot(ACTmath~semester,data=score)

# Trekker ut score for hoest 
score.hoest = score$ACTmath[score$semester==1]

# Ser paa oppsummmering og plot for hoestdataene:
summary(score.hoest)
qqnorm(score.hoest)
qqline(score.hoest)

# Trekker ut score for vaar
score.vaar = score$ACTmath[score$semester==2]

# Ser paa oppsummmering og plott for vaardataene:
summary(score.vaar)
qqnorm(score.vaar)
qqline(score.vaar)

# Antar at variansene i gruppe begge er like, sigma1=sigma2=4
# Får da 95% konfidensintervall ved følgende kommandoer

m = length(score.hoest)
n = length(score.vaar)
m1 = mean(score.hoest)
m2 = mean(score.vaar)
sigma1=sigma2=4

low = m1-m2-1.96*sqrt(sigma1^2/m+sigma2^2/n)
up = m1-m2+1.96*sqrt(sigma1^2/m+sigma2^2/n)
c(low,up)


# Lager så 95% konfidensintervall for forventet forskjell paa hoest og vaar basert på
# estimerte standardavvik
m = length(score.hoest)
n = length(score.vaar)
m1 = mean(score.hoest)
m2 = mean(score.vaar)
s1.2 = var(score.hoest)
s2.2 = var(score.vaar)
sqrt(c(s1.2,s2.2))
low = m1-m2-1.96*sqrt(s1.2/m+s2.2/n)
up = m1-m2+1.96*sqrt(s1.2/m+s2.2/n)
c(low,up)


## Vektoekning for rotter (igjen)

# Leser data
rotter = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/rotter.vekt.txt",header=TRUE)

# Boxplot
boxplot(vektok~beh,data=rotter)

# Trekker ut vektoekning for ozon-gruppa 
vekt.ozon = rotter$vektok[rotter$beh==1]

# Ser paa oppsummmering og plott for ozon-gruppa:
summary(vekt.ozon)
qqnorm(vekt.ozon)
qqline(vekt.ozon)

# Trekker ut vektoekning for kontrollgruppa 
vekt.kont = rotter$vektok[rotter$beh==2]

# Ser paa oppsummmering og plot for kontrollgruppa:
summary(vekt.kont)
qqnorm(vekt.kont)
qqline(vekt.kont)

# Lager 95% konfidensintervall for forventet forskjell paa ozon 
# og kontroll under forutsetning av lik varians i de to gruppene
# (Vi kan ogsaa bruke kommandoen "t.test"; se nedenfor)
m = length(vekt.ozon)
n = length(vekt.kont)
c(m,n)
m1 = mean(vekt.ozon)
m2 = mean(vekt.kont)
c(m1,m2)
s1.2 = var(vekt.ozon)
s2.2 = var(vekt.kont)
sqrt(c(s1.2,s2.2))
sp.2 = ((m-1)*s1.2+(n-1)*s2.2)/(m+n-2)
low.p = m1-m2-qt(0.975,m+n-2)*sqrt(sp.2*(1/m+1/n))
up.p = m1-m2+qt(0.975,m+n-2)*sqrt(sp.2*(1/m+1/n))
c(low.p,up.p)
t.test(vekt.ozon,vekt.kont,var.equal=TRUE)$conf.int

# Lager 95% konfidensintervall for forventet forskjell paa ozon 
# og kontroll uten aa forutsette lik varians i de to gruppene
# (Vi kan ogsaa bruke kommandoen "t.test"; se nedenfor)
df = (s1.2/m+s2.2/n)^2/((s1.2/m)^2/(m-1)+(s2.2/n)^2/(n-1))
low.s = m1-m2-qt(0.975,df)*sqrt(s1.2/m+s2.2/n)
up.s = m1-m2+qt(0.975,df)*sqrt(s1.2/m+s2.2/n)
c(low.s,up.s)
t.test(vekt.ozon,vekt.kont,var.equal=FALSE)$conf.int

## Eks. 10.8

# Leser dataene
exam = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp10-04.txt",header=TRUE)
head(exam)
summary(exam)
table(exam$GROUP)

# Boxplot 
boxplot(FINAL~GROUP,data=exam)

# Trekker ut score for tradisjonell undervisning (kontroll) 
exam.kontroll = exam$FINAL[exam$GROUP == "CONTROL"]

# Ser paa oppsummmering og plot for kontroll dataene:
summary(exam.kontroll)
hist(exam.kontroll)
qqnorm(exam.kontroll)
qqline(exam.kontroll)

# Trekker ut score for de som fikk nye undervisningsmetoder (eksperimentell) 
exam.exper = exam$FINAL[exam$GROUP=="EXPER"]

# Ser paa oppsummmering og plot for eksperimentelle data:
summary(exam.exper)
hist(exam.exper)
qqnorm(exam.exper)
qqline(exam.exper)

# Tester H_0: mu.kont >= mu.exper mot H_a: my.kont < my.exper
m = length(exam.kontroll)
n = length(exam.exper)
c(m,n)
m1 = mean(exam.kontroll)
m2 = mean(exam.exper)
s1.2 = var(exam.kontroll)
s2.2 = var(exam.exper)
c(m1,m2)
c(s1.2,s2.2)
z = (m1-m2)/sqrt(s1.2/m+s2.2/n)
p.value = pnorm(z)
c(z,p.value)

t.test(exam.exper,exam.kontroll,var.equal=T,alternative="greater")
t.test(exam.exper,exam.kontroll,var.equal=F,alternative="greater")

t.test(exam.exper,exam.kontroll,var.equal=T)
t.test(exam.exper,exam.kontroll,var.equal=F)

# Ved signifikansnivC% 5% forkaster vi H_0

# Vi proever aa fjerne studentene med de daarligste resultatene fra begge grupper
exam.kontroll.2 = exam.kontroll[exam.kontroll >= 10]
exam.exper.2 = exam.exper[exam.exper >= 10]
m = length(exam.kontroll.2)
n = length(exam.exper.2)
c(m,n)
m1 = mean(exam.kontroll.2)
m2 = mean(exam.exper.2)
s1.2 = var(exam.kontroll.2)
s2.2 = var(exam.exper.2)
c(m1,m2)
c(s1.2,s2.2)
z = (m1-m2)/sqrt(s1.2/m+s2.2/n)
p.value = pnorm(z)
c(z,p.value)

# NC% forkaster vi ikke H_0 ved signifikansnivaa 5%. Indikerer at er det for de
# svakeste studentene at de nye undervisningsmetodene gir utslag

## Vektoekning for rotter (forts.)

# Leser data og trekker ut de to gruppene
rotter = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/rotter.vekt.txt",header=TRUE)
vekt.ozon = rotter$vektok[rotter$beh==1]
vekt.kont = rotter$vektok[rotter$beh==2]

# T-test under forutsetning av lik varians
t.test(vekt.ozon,vekt.kont,var.equal=TRUE)

# T-test uten aa forutsette lik varians
t.test(vekt.ozon,vekt.kont)

## Teststyrke, feil av type II, utvalgsstoerrelse

# Vi ser paa testing av H_0: mu_1 = mu_2 mot H_a: mu_1 != mu_2 
# naar m = n = 20 og sigma = 15 

# Vi finner styrken til testen naar mu_1-mu_2 = delta
power.t.test(n=20,delta=10,sd=15,power=NULL,type="two.sample",alternative="two.sided")

# Vi finner m = n slik at testen faar styrke 90% naar delta = 10
power.t.test(n=NULL,delta=10,sd=15,power=0.90,type="two.sample",alternative="two.sided")





#----------------------------------------------------------------
# Parede data, Eks. 10.11
tidfør=c(81,87,86,82,90,86,96,73,74,75,72,80,66,72,56,82)
tidetter=c(78,91,78,78,84,67,92,70,58,62,70,58,66,60,65,73)
boxplot(tidfør,tidetter)

plot(tidfør,tidetter)
cor(tidfør,tidetter)

diff=tidfør-tidetter
boxplot(diff)
summary(diff);sd(diff)

t.test(diff)



## Eks. 10.11: muskel- og skjelettlidelser i nakke og skulder
#Tid observert hvor armheving var under 30 grader
#To observasjoner per pasient, 18 mnd imellom

# Leser dataene
arm.elevation=read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp10-09.txt",header=T)
head(arm.elevation)
summary(arm.elevation)

plot(arm.elevation$Before,arm.elevation$After)
cor(arm.elevation$Before,arm.elevation$After)

# Normalfordelingplott for differansene
qqnorm(arm.elevation$Diff)
qqline(arm.elevation$Diff)

# Test og konfidensinterval
t.test(arm.elevation$Diff,alternative="two.sided")


#------------------------------------------------
## Meningsmaaling fra Norstat 17. august 2021 og 17. august 2022.
## I 2022 svarte 41 av 691 at de ville stemme pC%
## SP hvis det var stortingsvalg. I 2021 svarte 117 av 677 
## at de vile stemme pC% SP. Har andelen som ville stemt SP endret
## seg signifikant?

# La p1 vaere SPs oppslutning i august 2021 og la p2
# vaere oppslutningen i august 2022. Vi vil lage et 
# konfidensintervall for p1-p2

m = 677
n = 691
x = 117
y = 41
p1.hat = x/m
p2.hat = y/n
c(p1.hat,p2.hat)
se = sqrt(p1.hat*(1-p1.hat)/m+p2.hat*(1-p2.hat)/n)
low = p1.hat-p2.hat-1.96*se
up = p1.hat-p2.hat+1.96*se
c(low,up)

# Vi vil ogsaa teste av H_0: p1 = p2 mot H_a: p1 != p2
# Beregner z-obsertator og P-verdi
p.hat = (x+y)/(m+n)
z = (p1.hat-p2.hat)/sqrt(p.hat*(1-p.hat)*(1/m+1/n))
p.value = 2*(1-pnorm(abs(z)))
c(z,p.value)

# Vi kan ogsaa bruke kommandoen "prop.test":
prop.test(c(x,y),c(m,n),alternative="two.sided",correct=F)
# prop.test angir test-observatir X-square=42.009=Z^2
sqrt(43.12)

# Med korreksjon for kontinuitet faar vi 
prop.test(c(x,y),c(m,n),alternative="two.sided")


## Eks. 10.13 Salk polio vaksine eksperiment 1954  

# Planlegging av eksperiment
# Vi lar p.1 og p.2 vaere sannsynlighetene for at et barn faar polio
# i henholdsvis kontroll- og behandlingsgruppa. Vi vil teste H_0: p.1 <= p.2 
# mot H_a: p.1 > p.2 med alpha = 0.05. Vi vil finne n (= m) slik at 
# sannsynligheten for Type II-feil er hoeyst beta naar p.1 = 0.0003 og p.2 = 0.00015.
p.1 = 0.0003
p.2 = 0.00015
alpha = 0.05
beta = 0.1
n = (qnorm(1-alpha)*sqrt((p.1+p.2)*(1-(p.1+p.2)/2))+qnorm(1-beta)*sqrt(p.1*(1-p.1)+p.2*(1-p.2)))^2/(p.1-p.2)^2
n

# For framtiden bruk software, f.eks.
power.prop.test(p1=p.1,p2=p.2,sig.level=alpha,power=1-beta,alternative="one.sided")


# Resultater fra det faktiske forsoeket
m = 201229
x = 110
n = 200745
y = 33
p.1.hat = x/m
p.2.hat = y/n
c(p.1.hat,p.2.hat)

p.hat = (x+y)/(m+n)
z.obs = (p.1.hat-p.2.hat)/sqrt(p.hat*(1-p.hat)*(1/m+1/n))
p.val = 1-pnorm(z.obs)
c(z.obs,p.val)

prop.test(x=c(x,y),n=c(m,n),correct=FALSE,alternative="greater")

# Resultatene gir overveldende stoette for aa forkaste H_0, og dermed
# for at Salks poliovaksine virker.



#_____________________________________________________________________________________________
# Leser data og trekker ut de to gruppene
rotter = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/rotter.vekt.txt",header=TRUE)
vekt.ozon = rotter$vektok[rotter$beh==1]
vekt.kont = rotter$vektok[rotter$beh==2]

# Vi vil lage et 95% konfidensintervall for forholdet 
# mellom variansen i ozon gruppen og kontrollgruppen

m = length(vekt.ozon)
n = length(vekt.kont)
s1.2 = var(vekt.ozon)
s2.2 = var(vekt.kont)
c(s1.2,s2.2)
alpha = 0.05
low = (s1.2/s2.2)*(1/qf(1-alpha/2,m-1,n-1))
up = (s1.2/s2.2)*(1/qf(alpha/2,m-1,n-1))
c(low,up)

# Vi tester nullhypotesen om like varianser
f = s1.2/s2.2
p.value = 2*(1-pf(f,m-1,n-1))
# Vi regner P-verdien paa denne maaten siden f > 1
# Hvis vi hadde observert f < 1, ville P-verdien vaert 2*pf(f,m-1,n-1)
c(f,p.value)

# Vi kan bruke kommandoen "var.test" til aa teste
# om variansene er like og beregne konfidensintervall
var.test(vekt.ozon,vekt.kont,alternative="two.sided")