################## ESIMERKKI 5.1 #################################

library(alr3)
data(longshoots)
attach(longshoots)
names(longshoots)

plot(Day, ybar, xlab="Päiviä kasvukauden alusta", 
ylab="Runkoyksiköiden keskiarvo")

lm.n <- lm(ybar ~ Day, weights=n)
abline(coef(lm.n))

summary(lm.n)

plot(fitted(lm.n), resid(lm.n), xlab="Sovite", ylab="Jäännös")

confint(lm.n)

################## ESIMERKKI 5.2 ###################################

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

attach(alko)
hinta <- ts(hinta, start=vuosi[1])
kulutus <- ts(kulutus, start=vuosi[1])
tulot <- ts(tulot, start=vuosi[1])
aika <- ts(vuosi-vuosi[1]+1, start=vuosi[1])
d69 <- ts(as.numeric(vuosi >= 1969), start=vuosi[1])

plot.ts(cbind(log(kulutus), log(hinta), log(tulot)), main="")

plot.ts(cbind(kulutus,hinta,tulot))

lm.alko <- lm(log(kulutus) ~ log(hinta) + log(tulot) + aika + d69)

round(rbind(coef(lm.alko),summary(lm.alko)$coef[,2]),4)

res <- ts(resid(lm.alko), start=vuosi[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.alko, test="Autocorrelation")

arima.alko <- arima(log(kulutus), order=c(2,0,0), 
xreg=cbind(log(hinta), log(tulot), aika, d69))

arima.alko

res.ar <- resid(arima.alko)
acf(res.ar)

 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,3,0)) 

 plot(res[-n], res[-1], xlab="", ylab="")


 mtext("Edellinen", font=1, side=1, line=2.5) 
 mtext("Jälkimäinen", font=1, side=2, line=2.5)
 mtext("Tavallinen regressio", font=1, side=3, line=1.5)

 par(mar=c(3,4,3,0))

 plot(res.ar[-n], res.ar[-1], xlab="", ylab="")


 mtext("Edellinen", font=1, side=1, line=2.5) 
 mtext("Jälkimäinen", font=1, side=2, line=2.5)
  mtext("Autoregressiivinen malli", font=1, side=3, line=1.5)

 par(oldpar)


plot(res.ar[-n], res.ar[-1], xlab="Edellinen", ylab="Jälkimmäinen")

#################### ESIMERKKI 5.3 #################################

plantgrowth <- read.table("http://users.jyu.fi/~junyblom/plantgrowth.dat",
header=T)
attach(plantgrowth)
names(plantgrowth)
# Käsittelyyn 2 liittyvät havainnot ruukuittain
cbind(growth[treat==2 & pot ==1],growth[treat==2 & pot ==2],
growth[treat==2 & pot ==3])

# Ruukkukeskiarvot
means <- tapply(growth, list(pot, treat), mean)
means

# Keskiarvot vektoreiksi, niihin liittyvät käsittelyt ja analyysi

means <- c(means)
treat.m <- factor(rep(1:6,each=3))
fixed.lm <- lm(means ~ treat.m)
summary(fixed.lm)

light.m <- factor(rep(rep(c(8,12,16), each=3), 2))
temp.m <- factor(rep(c("low","high"), each=9))

fixed2.lm <- lm(means ~ light.m + temp.m)
summary(fixed2.lm)

confint(fixed2.lm)

k <- 6; m <- 3; r <- 4
SSU <- sum(resid(fixed.lm)^2)

pot <- factor(pot)
treat <- factor(treat)

potmeans.lm <- lm(growth ~ pot*treat)
SSE <- sum(resid(potmeans.lm)^2)

# Estimaatit
sigma2.eps <- SSE/(k*m*(r-1))
sigma2.b <- (SSU/(k*(m-1))) - sigma2.eps/r
rho <- sigma2.b/sigma2.eps
c(sigma2.eps , sigma2.b, rho)

# Luottamusväli rho:lle ja testin p-arvo hypoteesille sigma^2_b = 0

up <- qf(.975, df1=k*(m-1), df2=k*m*(r-1))
lo <- qf(.025, df1=k*(m-1), df2=k*m*(r-1))
F <- ((r*SSU)/(k*(m-1))) / (SSE/(k*m*(r-1)))
1 - pf(F, df1=k*(m-1), df2=k*m*(r-1))
c((F/up -1)/r, (F/lo -1)/r)

# Analyysi lme() -funktiolla
library(nlme)

pot2 <- 10*as.numeric(pot) + as.numeric(treat)
pot2 <- factor(pot2)

growth.lme <- lme(growth ~ treat, random = ~ 1 | pot2)
summary(growth.lme)

pot.eff <- random.effects(growth.lme)[[1]]
qqnorm(pot.eff); qqline(pot.eff)

light <- factor(light)
growth.flme <- lme(growth ~ light + temp, random = ~ 1 | pot2)
summary(growth.flme)