## Teht 1
library(alr3)
data(allshoots)
attach(allshoots)
malli1<-lm(ybar~Day+Type,weights=n)
malli2<-lm(ybar~Day+Type+Type*Day,weights=n)
summary(malli2)
# Kahden ryhmän kulmakertoimet poikkeavat toisistaan
# Lyhyillä kulmakerroin = 0.19
# Pitkillä 0.22
# Kasvukauden alussa lyhyiden versojen runkoyksiköiden kappalemäärä on 9.5
# pitkillä 10.0
# Ero ei merkitsevä p = 0.19
plot(Day[Type==0], ybar[Type==0])
points(Day[Type==1], ybar[Type==1],col=2)
plot(malli2$fitted,malli2$resid)
plot(Day,malli2$resid)
plot(Type,malli2$resid)
boxplot(malli2$resid~Type)

qqnorm(malli2$resid)
qqline(malli2$resid)
# Lineaarisuus ja homoskedastisuus näyttävät toteutuvan tässä mallissa,
# selitysaste korkea
# Erilaiset Box-Cox - muunnokset vain huonontaisivat tilannetta



### TEHTÄVÄ 2 ####


earnings <-read.table("http://users.jyu.fi/~junyblom/earnings.dat", header=TRUE)
earnings$logearn <- log(earnings$earn)
earnings$earn14 <- (earnings$earn)^0.25
earnings$earn12 <- (earnings$earn)^0.5
earnings$earn34 <- (earnings$earn)^0.75 
attach(earnings)


# a)

malli1 <- lm(log(earn) ~ height+age+sex)  # lambda = 0 (logaritminen)
malli2 <- lm(earn14 ~ height+age+sex)     # lambda = 1/4
malli3 <- lm(earn12 ~ height+age+sex)	   # lambda = 1/2 
malli4 <- lm(earn34 ~ height+age+sex)	   # lambda = 3/4
malli5 <- lm(earn ~ height+age+sex)	   # lambda = 1

par(mfrow=c(3,2))
qqnorm(malli1$resid, main="Lambda = 0 (logaritminen)"); qqline(malli1$resid)
qqnorm(malli2$resid, main="Lambda = 1/4"); qqline(malli2$resid)
qqnorm(malli3$resid, main="Lambda = 1/2"); qqline(malli3$resid)
qqnorm(malli4$resid, main="Lambda = 3/4"); qqline(malli4$resid)
qqnorm(malli5$resid, main="Lambda = 1"); qqline(malli5$resid)

# qqnorm-plot on paras ("pisteet ovat lähimpänä suoraa"), kun lambda = 1/4


# b)

# Nyt lambda = 1/4
plot(malli2$fitted, malli2$resid, ylab="Jäännös", xlab="Sovite")
plot(height, malli2$resid, ylab="Jäännös", xlab="Pituus")
plot(age, malli2$resid, ylab="Jäännös", xlab="Ikä")

boxplot(malli2$resid~sex)
plot(malli2$fitted, malli2$resid, ylab="Jäännös", xlab="Sovite")
points(malli2$fitted[sex==1], malli2$resid[sex==1], col=3)
# Jäännös-sovite-kuviossa kaksi ryhmää (sukupuolen mukaan)
# Jäänös-ikä-kuviossa epälineaarisuutta "käänteinen-U" 

#####-------------------------------------------------------------
teht3

malli0<-lm(earn^.25~height + age +sex)
qqnorm(malli0$resi)
plot(malli0$fit,malli0$resi)
sex<-factor(sex)
a)

malli1<-lm(earn^.25~height + age +I(height^2) + I(age^2) +sex)
qqnorm(malli1$resi)
qqline(malli1$resi)
plot(malli1$fit,malli1$resi)
summary(malli1)
#iän neliö merkitsevä, pituuden ei.

malli2<-lm(earn^.25~height + age  + I(age^2) +sex)
summary(malli2)
qqnorm(malli2$resi)
plot(malli2$fit,malli2$resi)
#ei muutosta edelliseen

b)
malli3<-lm(earn^.25~height + age + I(age^2)+sex + 
sex*height+sex*age+sex*I(age^2))
summary(malli3)
#pituuden ja sukupuolen interaktio ei ole merkitsevä, 
#mutta iän- ja iän neliön interaktiot sukupuolen kanssa ovat.

c)
malli4<-lm(earn^.25~height + age + I(age^2)+sex +sex*age+sex*I(age^2))
summary(malli4)
# Miesten malli iän osalta vakio + 0.40*age - 0.0039*age^2
# Naisten malli iän osalta vakio + (0.40-0.19)*age -
# (0.0039 - 0.0018)*age^2
# ikä, missä maksimiansiot
b <- coef(malli3)
# Miehet
 -b["age"]/(2*b["I(age^2)"])
50.6
# Naiset
 -(b["age"]+b["age:sex2"])/(2*(b["I(age^2)"]+b["I(age^2):sex2"]))
49.9

### KERTOIMIEN TULKINNAT HANKALIA!!!!


## TEHT 4

salary <- read.table("http://users.jyu.fi/~junyblom/salary.dat",
header=T)
names(salary)
attach(salary)
Rank <- factor(Rank)
Degree <- factor(Degree)
Sex <- factor(Sex)

palkka_kaikki <- lm(Salary ~ Sex + Rank + Year + Degree + YSdeg)
summary(palkka_kaikki)
## 

palkka1 <- lm(Salary ~ Rank + Year + YSdeg + Degree)
summary(palkka1)

plot(Sex, Salary, xlab="Sukupuoli", ylab="Palkka")
## Aineistossa yksi suuri palkkainen nais professori

palkka0 <- lm(Salary ~ Rank + Year)
summary(palkka0)

anova(palkka0,palkka_kaikki, test="F")

## Palkka0 Riittävä malli, virka-aseman kasvaessa , palkka kasvaa. 
## Samassa virka-asemassa ollessa palkka kasvaa keskimäärin 375 dollaria vuodessa.

# Tehtävä 5

spirit <- read.table("http://users.jyu.fi/~junyblom/spirit.dat",
header=T)
attach(spirit)

spirit10 <- 10^spirit
logspirit <- log(spirit10)
time <- 1:69

spirit.uusi <- data.frame(logspirit,time)

consumption <- spirit.uusi$consumption
income <- spirit.uusi$income
price <- spirit.uusi$price
spirit.uusi

price <- ts(price, start=time[1])
consumption <- ts(consumption, start=time[1])
income <- ts(income, start=time[1])
time <- ts(time-time[1]+1, start=time[1])

plot.ts(cbind(consumption, price, income), main="")

lm.spirit <- lm(consumption ~ price + income + time)
summary(lm.spirit)

res <- ts(resid(lm.spirit), start=time[1])
n <- length(res)
plot(res[-n],res[-1], xlab="Edellinen", ylab="Jälkimmäinen")

sum(res[-n]*res[-1])/sum(res^2)

source("http://users.jyu.fi/~junyblom/diagtest.R")
diagtest(lm.spirit, test="Autocorrelation")