# Testimisest:

x=rep(1:4, c(10, 20, 25, 28))
table(x)
y=c(0, 1, 5, -1)[x]+rt(length(x), df=1)

mudel=lm(y~as.factor(x))
anova(mudel)

kruskal.test(y~factor(x))

y2=rt(length(x), df=1)

mudel=lm(y2~as.factor(x))
anova(mudel)
kruskal.test(y2~factor(x))

kordi=10000
tul=rep(NA, kordi)
tul2=rep(NA, kordi)
for (i in 1:kordi){
  y=c(0, 1, 2.5, -1)[x]+rt(length(x), df=1)
  mudel=lm(y~as.factor(x))
  tul[i]=anova(mudel)[1,5]
  tul2[i]=kruskal.test(y~factor(x))$p.value
}



# ----------------------------------------------------------------
# Tihedusfunktsiooni hindamisest ---------------------------------
# ----------------------------------------------------------------

# Algandmed
x=c(1, 10, 12, 20, 40, 45, 46, 55)

# Hindame tihedusfunktsiooni kasutades tummafunktsioonina K()
# normaaljaotuse tihedust dispersiooniga 2:
plot(density(x, bw=2))
# Lisame joonise allserva tegelikud vaatlused
rug(x, lwd=3)

# Vaatame ühe vaatluse "kontributsiooni":
abix=seq(-10, 70, length=200)
lines(abix, 1/8*dnorm(abix-x[1], sd=2), col=2, lwd=2)

# Muudame kasutatava tuumafunktsiooni (tihedusfunktsiooni K(.)) standardhälvet
abi=density(x, bw=10)
lines(abi$x, abi$y, lwd=2, col="blue")

# Vaatame ühe vaatluse "kontributsiooni":
lines(abix, 1/8*dnorm(abix-x[1], sd=10), col="blue", lwd=2)


# Arvutame tiheduse hinnangu käsitsi:
lines(abix, 1/8*(
  dnorm(abix-x[1], sd=10)+dnorm(abix-x[2], sd=10)+
  dnorm(abix-x[3], sd=10)+dnorm(abix-x[4], sd=10)+
  dnorm(abix-x[5], sd=10)+dnorm(abix-x[6], sd=10)+
  dnorm(abix-x[7], sd=10)+dnorm(abix-x[8], sd=10)
    ), col=2, lwd=3, lty=2)


# Märkus: võib kasutada ka järgmist varianti (suure valimi puhul mugavam):

# valimi suurus
n=length(x)

# kui paljudes kohtades soovime arvutada tihedusfunktsiooni väärtust
narvutusi=length(abix)

tihedus=rep(0, narvutusi)
for (i in 1:n){
  tihedus=tihedus+dnorm(abix-x[i], sd=10)
}
tihedus=tihedus/n

plot(abix, tihedus, type="l")
rug(x)




# Oletame, et lisandub veel üks vaatlus (x=70):
x=c(x,70)

# Joonista uue, täiendatud andmestiku jaoks tihedusfunktsiooni hinnang.
# Kasuta tuumafunktsioonina t-jaotust (vabadusastmete arv=2, teisenda nii, et
# tuumafunktsiooni standardhälve oleks 10).
# t-jaotuse tihedusfunktsioon on R'is leitav käsuga dt(x, df=df).
# Võrdle saadud tulemust normaaljaotust tuumafunktsioonina kasutava hinnanguga 
# (kasuta mõlemal juhul sama tuumafunktsiooni hajuvust)


# -----------------------------------------------------------------
# Mitteparameetriline regressioon ---------------------------------
# -----------------------------------------------------------------

# Algandmed
data(cars); attach(cars)
plot(speed, dist)

# Leiame kaks mudelit - lokaalseid polünoome kasutades ja tavalise lin. reg:
mpreg = loess(dist ~ speed)
linreg = lm(dist~speed)

# Arvutame mõlema mudeli järgi prognoosid:
x= seq(0, 30, length=100)
ylin=predict(linreg, data.frame(speed=x))
ymp=predict(mpreg, data.frame(speed=x))

# Ja kanname leitud prognoosid regressioonijoontena joonisele:
lines(x, ylin)
lines(x, ymp, col="red", lwd=2)

# Lisame loess-joonele ligikaudse usaldusintervalli
mpse=predict(mpreg, data.frame(speed=x), se=TRUE)$se.fit
lines(x, ymp+1.96*mpse, lty=2, col="red")
lines(x, ymp-1.96*mpse, lty=2, col="red")

# Muudame hinnangut - lubame vaid 1/3 andmetest mõjutada regressioonijoont:
mpreg2 = loess(dist ~ speed, span=1/3, degree=1)
ymp2=predict(mpreg2, data.frame(speed=x))
lines(x, ymp2, col="blue", lwd=2)


# Proovime ise - apelsiinipuude peal:

# Algandmed
data(Orange)
attach(Orange)
plot(age, circumference)

# Leiame regressioonijoone:
n=100
prognoos=rep(NA, n)
x=seq(0, 2000, length=n)
for (i in 1:n){
  mudel=lm(circumference~age, weight=dnorm(age-x[i], sd=175))
  prognoos[i]=predict(mudel, data.frame(age=x[i]))
}
lines(x, prognoos, lwd=2)


# Põhimõtte illustratsioon - regressioonijoone leidmine punktis age=750:

abline(v=750, lty=2)
mudel=lm(circumference~age, weight=dnorm(age-750, sd=175))
y=predict(mudel, data.frame(age=x))
lines(x,y, col=2)
kaal=dnorm(age-750, sd=175)
kaal
# Standardiseerime kaalud vahemikku 0...1:
ukaal=(kaal-min(kaal))/(max(kaal)-min(kaal))
ukaal
# Näitame välja, millised punktid kui tugevalt mõjutavad 
#	regressioonijoone hinnangut punktis age=750
points(age, circumference, col=rgb(1-ukaal, 1-ukaal, 1-ukaal))
# Lisame kõik vaatlused imepisikeste täpikestena
points(age, circumference, pch=".")




# Muudame veidi programmi, arvutame kaalud tsipa teisiti - 
#	nimelt muudame akna laiust vastavalt teiste "naaber"vaatluste kaugusest.
# 	Muudame ka kasutatavat tuumafunktsiooni - teeme ta selliseks, et
#	kaugematel vaatlustel puuduks absoluutselt igasugune mõju...

plot(age, circumference)

# Kui suurt osa "naabervaatluseid" üldse kasutame
osak=0.5 
prognoos=rep(NA, n)
x=seq(0, 2000, length=n)
for (i in 1:n){
  kaugused=abs(x[i]-age)
  kvantiil=quantile(kaugused, osak)
  z=kaugused/(kvantiil)
  kaal=rep(0, length(z))
  kaal[z<1]=(1-abs(z[z<1])**3)**3
  mudel=lm(circumference~age, weight=kaal)
  prognoos[i]=predict(mudel, data.frame(age=x[i]))
}
lines(x, prognoos, lwd=2)

# Lisame joonisele punase värviga loess-kõvera 
mpreg = loess(circumference ~ age, span=0.5, degree=1, control=loess.control(surface = c("direct"), statistics=c("exact")))
ymp=predict(mpreg, data.frame(age=x))
lines(x, ymp, col="red", lwd=2)
# Mida näed?

# Uurime kasutatud kaale kohas age=750:
abline(v=750, lty=2)
kaugused=abs(750-age)
kvantiil=quantile(kaugused, osak)
z=kaugused/(kvantiil)
kaal=rep(0, length(z))
kaal[z<1]=(1-abs(z[z<1])**3)**3
kaal
points(age, circumference, col=gray(1-ukaal))

# Joonistame välja tuumafunktsiooni (x=750 jaoks):
vanus=seq(0, 2000, length=200)
kaugused=abs(vanus-750)
z=kaugused/(kvantiil)
kaal=rep(0, length(z))
kaal[z<1]=(1-abs(z[z<1])**3)**3

plot(vanus, kaal, type="l", main="K(.) kohas age=750")
abline(v=750, lty=2)



# -------------------------------------------------
# Muu: gruppide võrdlemine. Permutatsioonitest vs Kruskal-Wallis
# -------------------------------------------------

# Abifunktsioon - 
# Leiab F-statistiku väärtuse
Fstat = function(x, grupp)
{
	abi = table(grupp)
	SST = sum(abi * by(x, grupp, mean)^2) -length(x)*(mean(x))^2
	SST / (length(abi)-1) / ( (sum((x-mean(x))^2) - SST)/(length(x)-length(abi)))
}


# Algandmed:
x = rexp(9)
grps = c(1,1,1,2,2,2,3,3,3)
# x[grps==1] = x[grps==1] + 2
by(x, grps, mean)
boxplot(split(x,as.factor(grps)))


Fstat(x,grps)
Fobs = Fstat(x,grps)

permapproxF = function(x,grps,n)
{
	tulem = rep(NA,n)
	for (i in 1:n) tulem[i] = Fstat(sample(x),grps)
      stat = Fstat(x,grps)
      pvalue = mean(tulem >= stat)
	return(list(tulem=tulem,stat=stat,pvalue=pvalue))
}
ajut=permapproxF(x,grps,1000)
ajut$pvalue
# Olulisustõenäosuse hinnang on kirjas "mean of x" all, usaldusintervall iseloomustab 
#	olulisustõenäosuse arvutamise täpsust
t.test(1*(ajut$tulem>=Fobs))

#Kruskal-Wallis test
kruskal.test(x,grps)

# Muuda programmi selliselt, et grupp 1-te kuuluvatel vaatlustel oleks x kahe võrra suurem.
# 	Kumb meetoditest annab väiksema p-value? 


# Katsetame pärisandmetel:

## Vähihaigete elulemus
survtime = read.table(url("http://www.ms.ut.ee/mart/mpara2009/survtime.txt"),
                     col.names=c("type","gender","age","days","control","posttreat1","posttreat2"))
x <- survtime[,4] 
grps <- survtime[,1]  #1=stomach, 2=bronchus, 3=colon, 4=rectum, 5=bladder, 6=kidney

by(x,grps, mean)
table(grps)
boxplot(x~grps)

# tavaline dispersioonanalüüs
m1=aov(x ~ factor(grps))
summary(m1)

# Jäägid päris koledad....
qqnorm(resid(m1)); qqline(resid(m1))



# Permutatsioonitest - teststatistik sama mis dispersioonanalüüsil
Fobs = Fstat(x,grps)
Fobs

n=10000
abi = permapproxF(x,grps,n)
abi$pvalue
t.test(abi$tulem >= Fobs)
binom.test(sum(abi$tulem >= Fobs), n)

# Kruskal-Wallis, kasutades hii-ruut lähendit
kruskal.test(x,grps)


# -------------------------------------------------
# Veel permutatsioonitestist
# -------------------------------------------------
# Permutatsioonitest. Negatiivne(?) näide - halvasti
#   permuteerides võime testida midagi muud kui ise arvame:

kordusi=1000

pvalue=rep(NA, kordusi)
for(j in 1:kordusi){
# Kaks valimit populatsioonidest, mille hajuvus on sama
x1=rnorm(4, sd=2)
x2=rnorm(6, mean=21, sd=2)
x=c(x1,x2)
# Teststatistik - hajuvuste suhe - kontrollib hajuvuste erinevust?
stat=sd(x1)/sd(x2)

kord=100
statX=rep(NA,kord)
for (i in 1:kord){
grupp=sample(rep(0:1, c(4,6)))
xx1=x[grupp==1]
xx2=x[grupp==0]
statX[i]=sd(xx1)/sd(xx2)
}

pvalue[j]=sum(stat>statX)/kord
}
t.test(1*(pvalue<0.05))