## Teht1
earnings <- read.table("http://users.jyu.fi/~junyblom/earnings.dat",
header=T)
earnings$sex<-as.factor(earnings$sex)
attach(earnings)

plot(age,earn)

malli1<-lm(log(earn)~sex+age+height+sex*age,data=earnings)

summary(malli1)
# Sama sukupuoli, sama ikä, pituudessa 1 cm ero => ansio 1,0 % ero
# Miehet, sama pituus, iässä vuoden ero => ansio 1,3% ero
# Naiset, sama pituus, iässä vuoden ero => ansio 0,4% ero

k<-round(exp(coef(malli1)),5)
k

k[3]*k[5]

malli10<-lm(log(earn)~sex+age+height+sex*age+sex*height,data=earnings)
anova(malli10,malli1)
summary(malli10)

## Teht2

plot(height,log(earn))
plot(log(age),log(earn))
# Riippuvuus iän kanssa on ehkä hieman epälineaarista

plot(malli1$fit,malli1$resid)
# Heteroskaedastisuutta hsvaittavissa

qqnorm(malli1$resid)
qqline(malli1$resid)
# Jonkin verran epänormaaliutta, paljon pieniä arvoja

# Kuvassa näkyy ehkä joitakin poikkeuksellisen matalia arvoja
which.min(rstudent(malli1))

malli12<-lm(log(earn)~sex+age+height+sex*age,
data=earnings[rstudent(malli1)!= min(rstudent(malli1)),])
summary(malli12)
# Pieni vaikutus kertoimiin

max(hatvalues(malli1))
plot(hatvalues(malli1))
malli13<-lm(log(earn)~sex+age+height+sex*age,
data=earnings[-(hatvalues(malli1)==max(hatvalues(malli1)))])

summary(malli13)
# Ei suurta vaikutusta


#demo7 tehtävä 3
mesquite <- read.table("http://users.jyu.fi/~junyblom/mesquite.dat", header=TRUE)
attach(mesquite)

can.area <- diam1*diam2

# Malli C
lm.c0 <- lm(log(leafwt) ~ log(can.area) + log(canht) + group)

lm.c<-lm(log(leafwt) ~ log(can.area) + log(canht) + group+log(canht)*log(can.area)*group)
anova(lm.c0,lm.c,test="F")

lm.c<-lm(log(leafwt) ~ log(can.area) + log(canht) + group+log(canht)*log(can.area))
anova(lm.c0,lm.c,test="F")

lm.c<-lm(log(leafwt) ~ log(can.area) + log(canht) + group+log(canht)*group)
anova(lm.c0,lm.c,test="F")

lm.c<-lm(log(leafwt) ~ log(can.area) + log(canht) + group+log(can.area)*group)
anova(lm.c0,lm.c,test="F")

#Missään tapauksessa interaktio ei ole merkitsevää

####################TEHT 4 ##############################################

library(foreign)
kidiq = read.dta("http://www.stat.columbia.edu/~gelman/arm/examples/child.iq/kidiq.dta")
attach(kidiq)
names(kidiq)
mom_hs = factor(mom_hs)
mom_work = factor(mom_work)

#Malli1
malli1 = lm(kid_score ~ mom_hs + mom_iq + mom_work + mom_age + mom_hs*mom_iq + mom_iq*mom_work)
anova(lm(kid_score ~ mom_hs + mom_iq + mom_work + mom_age + mom_hs*mom_iq),malli1, test="F")
summary(malli1)
#mom_work*mom_iq pois

#Malli2
malli2 = lm(kid_score ~ mom_hs + mom_iq + mom_work + mom_age + mom_hs*mom_iq)
summary(malli2) 
#age pois

#Malli3
malli3 = lm(kid_score ~ mom_hs + mom_iq + mom_work + mom_hs*mom_iq)
summary(malli3) 
anova(lm(kid_score ~ mom_hs + mom_iq + mom_hs*mom_iq),malli3,test="F")
#mom_work pois

#Malli4
malli4 = lm(kid_score ~ mom_hs + mom_iq + mom_hs*mom_iq)
summary(malli4)

#verrataan äitejä, jotka eivät ole käyneet high schoolia ja äö 
#poikkeaa pisteellä, niin lasten testipistemäärät poikkeavat 0.97 pisteellä

#verrataan äitejä, jotka ovat käyneet hs:n ,niin lasten testipistemäärät
#poikkeavat 0.48 pistettä




################ TEHT 5 ########################################

malli = lm(kid_score ~ mom_hs + mom_iq + mom_hs*mom_iq)
plot(fitted(malli1),resid(malli1))
#lineaarisuus ja heteroskedastisuus OK
qqnorm(resid(malli1)); qqline(resid(malli1))
#ehkä pienet jäännökset hiukan liian pieniä,
#normaalisuus ~Ok

#Outliereita?
which.max(resid(malli))
malli2 = lm(kid_score ~ mom_hs + mom_iq + mom_hs*mom_iq,subset=(87 != 434))
malli
malli2

#Vipuarvoja?
plot(hatvalues(malli), xlab="Tapaus", ylab="Vipuarvo")
#alle 0.1 -> ei ongelma