andmed=read.csv("http://haldna.ee/r/Statistika/EMU_2023.csv",header=TRUE,sep=",", dec=".")
attach(andmed)

names(andmed)
andmed

summary(andmed)

boxplot(pikkus,ylab="Pikkus, cm",main="Kõik tudengid")
boxplot(kaal,ylab="Kaal, kg",main="Kõik tudengid")
mean(pikkus)
sd(pikkus)
hist(pikkus)

Pikkus=rnorm(n=300,mean=172,sd=9)

par(mfrow=c(1,2)) # kaks graafikut kõrvuti
hist(pikkus, main="Empiirilised ",ylab="tudengite arv")
hist(Pikkus, main="Teoreetilised",ylab="")

kaal2=rnorm(n=300,mean(kaal),sd(kaal))
hist(kaal, main="Empiirilised kaalud",ylab="tudengite arv")
hist(kaal2, main="Teoreetilised kaalud",ylab="")

table(mathinne)
prop.table(table(mathinne)) # suhtelised sagedused
prop.table(table(mathinne))*100 # suhtelised sagedused protsentides
round(prop.table(table(mathinne))*100,1)

T1=table(mathinne)
prop.table(T1)*100
round(prop.table(T1)*100,1)
barplot(T1)
pie(T1) # võimalikud graafikud

shapiro.test(pikkus) 

Kas naiste ja meeste pikkuste varieeruvused on sarnased ?
bartlett.test(pikkus~factor(sugu))
fligner.test(pikkus~factor(sugu))
var.test(pikkus~factor(sugu))
table(sugu)

sd(pikkus[sugu=="m"])
sd(pikkus[sugu=="n"])
names(andmed)

boxplot(pikkus~factor(sugu),col=FALSE)

boxplot(biolpunkte~factor(aeg),col=FALSE)

boxplot(kaal~factor(aeg),col=FALSE)

t.test(pikkus~factor(sugu))

t.test(pikkus)

# Teine praktikum

andmed=read.csv("http://haldna.ee/r/Statistika/EMU_2023.csv",header=TRUE,sep=",", dec=".")
attach(andmed)

names(andmed)

T1=table(sugu,aeg);addmargins(T1)

chisq.test(T1)

barplot(T1,beside=TRUE,legend=TRUE)
barplot(prop.table(T1,margin=2),legend=TRUE,ylab="Tudengite osakaal",xlim=c(0,4.5))

prop.test(c(38,42),n=c(50,150))

TT<-matrix(c(50,25,45,80),nrow=2,byrow=TRUE,dimnames=list(c("Suvila on","Suvilat pole"),c("Taime pole","Taim on")));TT

fisher.test(TT) # kas suhe tuli sama, mis käsitsi arvutades?

barplot(TT,legend=TRUE,xlim=c(0,3.5))

barplot(t(TT),legend=TRUE,xlim=c(0,3.5))
barplot(prop.table(t(TT),margin=2)*100,legend=TRUE,ylab="Taime olemasolu,%",xlim=c(0,3.5))
chisq.test(TT)

T2=table(sugu,aeg) ;T2
fisher.test(T2)
table(aeg)
summary(andmed)
table(tase)

t.test(biolpunkte~factor(tase))
t.test(pikkus~sugu)
names(andmed)
kor_andmed=andmed[,-c(1,4,10:13)]

maatriks=round(cor(kor_andmed, use="pairwise.complete.obs"),2)
maatriks

cor.test(biolpunkte,inglise,method="pearson")

cor.test(pikkus,kaal,method="pearson")
cor.test(pikkus,kaal,method="spearman")
cor.test(pikkus,kaal) # mis see on?


Iseseisev töö, kirjuta välja eelneva näite H0 ja H1.

plot(pikkus,kaal)

plot(kaal~pikkus)

# vaatame kolme erinevat joonist kõrvuti
par(mfrow=c(1,3))
plot(pikkus,kaal,main="plot(pikkus,kaal)")
plot(kaal,pikkus,main="plot(kaal,pikkus)")
plot(kaal~pikkus,main="plot(pikkus~kaal)")
#Lisame ühikud telgedele
plot(kaal~pikkus,main="Kaalu seos pikkusega",ylab="Kaal,kg",xlab="Pikkus, cm")
#muudame graafiku akent, kui vaja
pairs(kor_andmed)
#install.packages("corrplot")
library(corrplot) 
#Korrelatsioonimaatriksi graafiliseks esitamiseks
corrplot(maatriks)

# Lineaarsed mudelid

fosfor<-read.table("http://ph.emu.ee/~haldna/Statistika/CY_naide2.csv",
sep=",", dec=".", header=TRUE, na.strings=c("NA", ""))

names(fosfor);attach(fosfor)
summary(fosfor)
punkt=factor(punkt)
summary(punkt)
table(punkt)
table(PTOT)

# kontrolli, millised PTOT ja PO4P väärtused on üldse andmebaasis olemas,
# vaata ka jaotust!
hist(PTOT)
shapiro.test(PTOT)

hist(PO4P) # kas on n.j.?

cor.test(PTOT,PO4P) # kas lin.korrelatsioon on oluline
plot(PTOT~PO4P)

lin_mudel=lm(PTOT~PO4P)

summary(lin_mudel)

# arvutame prognoosid paberil käsitsi või kalkulaatori abil
# üldfosfor=....+....*fosfaat
# R abil

# kui PO4P on 25
 32.3714+1.7340*25

# Joonis, kõigepealt hajuvusgraafik uuesti
plot(PTOT~PO4P)
abline(lin_mudel)
#Prognooside leidmiseks teeme valiku argumenttunnuse (PO4P) väärtustest.
# Nendele leiab R on prognoosi, usaldusvahemiku ja prognoosi vahemiku.
# Lihtsa regressioonimudeli korral piisab kahest etteantud väärtusest,
# võtame näiteks minimaalse ja maksimaalse PO4P ja lisaks 10, kontrollimaks, kas arvutasime
#eespool õigesti.

uued=c(min(PO4P),10,max(PO4P))
uued
prognoos=predict(lin_mudel,data.frame(PO4P=uued))
prognoos #vaatame prognoose

us_piirid=predict(lin_mudel,data.frame(PO4P=uued),interval="conf")
us_piirid # vaatame prognoose koos 95% usalduspiiridega

pr_piirid=predict(lin_mudel,data.frame(PO4P=uued),interval="pred")
pr_piirid # vaatame prognoose koos 95% prognoosi vahemikega

# Joonis, kõigepealt hajuvusgraafik uuesti
plot(PTOT~PO4P)
abline(lin_mudel)
# NB!abline(lin_mudel) annab sama, kuid seda tehakse vaid ühe faktoriga #
#mudeli korral
lines(uued,prognoos)
lines(uued,us_piirid[,2],lty=2);lines(uued,us_piirid[,3],lty=2)
lines(uued,pr_piirid[,2],lty=3);lines(uued,pr_piirid[,3],lty=3)

####    Logaritmitud PTOT 

par(mfrow=c(1,2))
hist(PTOT);shapiro.test(PTOT)

hist(log(PTOT));shapiro.test(log(PTOT))

log_mudel=lm(log(PTOT)~PO4P);summary(log_mudel) # R2 on suurem

exp(3.472950+0.030607*10) #see on summary väljatrükist, piisav edasiseks
exp(3.472950)*exp(0.030607*10) # siin teisendame edasi

# prognoosid ja joonis algskaalas
uued2=c(0,10,20,40,50,60) # võtame rohkem X ehk PO4P väärtusi
prognoos2=exp(predict(log_mudel,data.frame(PO4P=uued2))) # tagasiteisendus
prognoos2
us_piirid2=exp(predict(log_mudel,data.frame(PO4P=uued2),interval="conf"))
pr_piirid2=exp(predict(log_mudel,data.frame(PO4P=uued2),interval="pred"))
plot(PTOT~PO4P)
lines(uued2,prognoos2,col=2) # siin ei kasuta enam abline
lines(uued2,us_piirid2[,2],lty=2);lines(uued2,us_piirid2[,3],lty=2)
lines(uued2,pr_piirid2[,2],lty=2);lines(uued2,pr_piirid2[,3],lty=2)

#Normaaljaotuse kontroll
hist(residuals(log_mudel))

shapiro.test(residuals(log_mudel))
shapiro.test(residuals(lin_mudel))

hist(CY);hist(log(CY))# kas on n.j.?
plot(log(CY)~d1)

mudel1<-lm(log(CY)~d1);summary(mudel1)
abline(mudel1)

hist(residuals(mudel1))
shapiro.test(residuals(mudel1))
plot(d1,residuals(mudel1), main=" Lineaarse mudeli jäägid")
lines(lowess(d1,residuals(mudel1)))

mudel2<-lm(log(CY)~ poly(d1, 2))
summary(mudel2)

mudel3=lm(log(CY)~poly(d1,3));summary(mudel3)
anova(mudel2,mudel3, test="Chisq")


#mudel3 pole parem

x=seq(1, 3.65, length=100)

y1=predict(mudel1, data.frame(d1=x), type="response")
y2=predict(mudel2, data.frame(d1=x), type="response")
#Kanname leitud prognoosid juba valmisolevale joonisele lines-käsu abil:
plot(d1, log(CY))
lines(x,y1, lwd=1, col=1)
lines(x,y2, lwd=1, col="red")
# teisendame CY prognoosid algskaalasse tagasi exp-funktsiooniga
plot(d1, CY)
lines(x,exp(y1), lwd=1, col=1)
lines(x,exp(y2), lwd=1, col="red")

max_mudel=lm(log(CY)~punkt+kuu+poly(d1,2)+NTOT+PTOT+CAE+MGE+O2+TEMPproov+PH_LF+ZB+SD) 
summary(max_mudel)

max2=step(max_mudel,direction="both",test="Chi")

summary(max2) # parim mudel
hist(residuals(max2))
shapiro.test(residuals(max2))

## logistiline regressioon

tigu<-read.csv("http://haldna.ee/r/Statistika/tigu.csv")
tigu
attach(tigu)
expcoef <- 
function(modelfit) {
OR <- exp(coef(modelfit))
cov.coef <- summary(modelfit)$cov.unscaled # covariance matrix
p <- summary(modelfit)$coef[,4]
SE.coef <- sqrt(diag(cov.coef))
CI95.lo = OR*exp(-1.96*SE.coef)
CI95.up = OR*exp(+1.96*SE.coef)
cbind(OR=round(OR,2),CI95.lo=round(CI95.lo,2),CI95.up=round(CI95.up,2),p=round(p,4))
 }

pyes=cbind(suri,elus); pyes
mud=glm(pyes~Liik+Moju+Niiskus+Temp,family=binomial)
summary(mud)

anova(mud,test="Chisq") # testime, millised faktorid on olulised
expcoef(mud)

# võtame aktiivseks faili fosfor

attach(fosfor)
table(dino1)
table(punkt)

mud_dino=glm(dino1~TEMPproov+FBM+punkt,family=binomial)
summary(mud_dino)

mud_dino2=glm(dino1~punkt,family=binomial)
summary(mud_dino2)
exp(1.2939)
exp(-1.2939)
expcoef(mud_dino2)
1/0.27

Dino_tabel=table(dino1,punkt)
addmargins(Dino_tabel)

fisher.test(Dino_tabel)

barplot(Dino_tabel,legend=TRUE)
barplot(prop.table(Dino_tabel,margin=2)*100,legend=TRUE)

# Kui palju väiksem on sanss (tõenäosus) leida "dino" punktis 16h võrreldes punktiga 11h
1/0.27

####   Dispersioonanalüüs   #####

data=read.csv("http://haldna.ee/r/Statistika/fosfor_disp.csv", sep=",", dec=".", header=TRUE) 

attach(data) 
data
interaction.plot(P,N,saagikus) 

anova_mud=aov(saagikus~factor(N)+factor(P)+(factor(N)*factor(P))) 

summary (anova_mud) 

TukeyHSD(anova_mud) 

model.tables(anova_mud,"means") 


lm_mud=lm(saagikus~factor(N)+factor(P)+(factor(N)*factor(P))) 

summary (lm_mud) 


Kalad=read.csv("http://haldna.ee/r/Statistika/kalasaagid.csv", sep=",", dec=".", header=TRUE, na.strings=c("NA", "")) 

attach(Kalad) 

names(Kalad);summary(Kalad)  # tutvun andmetega 

# mitu taset on faktortunnustel ja kui palju on igas mõõtmisi? 

 

table(liik);table(aasta);table(koht)  

 

# Vaatame faktorite võimalikku mõju püügikogusele (tunnus pyyk) graafiliselt. Kas mediaanid ja kvartiilid on erinevad?  

 

boxplot(pyyk~liik) 

boxplot(pyyk~aasta) 

boxplot(pyyk~liik*aasta) 

# Tee ise joonis faktori koht mõju vaatamiseks. 

 # Ühefaktoriline dispersioonanalüüs. Mudelis aasta faktor. 

 mud1a=aov(pyyk~factor(aasta));summary(mud1a)  # aasta mõju 

mud1l=aov(pyyk~liik);summary(mud1l)   # liikide erinevus 


# prognoosimise näide 

 

predict(mud1l, data.frame(liik=c("liik_1","liik_2","liik_3","liik_4")),interval="confidence") 

model.tables(mud1l,"means") 

# aritmeetilised keskmised liikide kaupa 

akesk=tapply(pyyk,liik,mean) # prognoosid langevad kokku  
akesk
sdh=tapply(pyyk,liik,sd)  # standardhälbed 

# Kahefaktoriline dispersioonanalüüs. Mudelis aasta ja liik. 

 mud2=aov(pyyk~factor(aasta)+liik);summary(mud2) 

interaction.plot(aasta,liik,pyyk) 
TukeyHSD(mud2)

 #### Kordamine
# kõigepealt salvestan kordamise andmed (hiirega parempoolne klõps->Save as->

tud=read.csv(file.choose(),sep=";",dec=",")

names(tud)
attach(tud)

#1. Kui palju on mehi, kui palju naisi?

  table(sugu)

#4.1 Mis on küsitletud tudengite matemaatika punktide (mate_punktid) keskmine väärtus (sobivad nii aritmeetiline keskmine kui ka mediaan)?  

#4.2. Kas  tunnus mate_punktid on aritmeetilist keskmist ja mediaani võrreldes normaaljaotusele lähedane (edaspidi n.j)?
summary(mate_punktid)

4.3. Kas tunnus mate_punktid on histogrammi vaadates n.j?
4.4. Kas tunnus mate_punktid on shapiro testi alusel n.j.?

hist(mate_punktid);shapiro.test(mate_punktid)

5. Leia valimi põhjal matemaatika punktide (mate_punktid) keskmise 95% usalduspiirid üldkogumi jaoks ( vihje t.test!)
t.test(mate_punktid)

5.2. Leia soo kaupa tudengite matemaatika punktide keskmised. ( vihje t.test!)

t.test(mate_punktid~sugu)

6.1. Kasuta eelmise hüpoteesi kontrollimiseks mitteparameetrilist meetodit (vihje: wilcoxon.test)
wilcox.test(mate_punktid~sugu)

6.2.  tee vastav joonis kasutades karp-vurrud diagrammi.
boxplot(mate_punktid~sugu)

7. Uurige kas matemaatika lõpueksami punktide keskmised on erinevad seoses sellega, kas tudeng käib trennis  või mitte (faktori treening). 
boxplot(mate_punktid~treening)

8.Millise kallakuga tudengeid on valimis kõige rohkem? Kasutage sagedustabelit, tulp- ja sektordiagrammi. (barplot, pie)
table(kallak)
prop.table(table(kallak))*100
pie(table(kallak))