Luku 1

################ ESIMERKKI 1.1 ###########################

read.table("http://users.jyu.fi/~junyblom/birthweight.dat",
 header=TRUE)

 attach(birthweight)

 plot(gestation[sex == 0], weight[sex ==0], xlab="Raskauden kesto",
 ylab="Syntymäpaino")

 out <- lm(weight~gestation, subset=(sex==0))
 coef(out)

 abline(coef(out))

################# ESIMERKKI 1.2 ############################

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

 pairs(fuel[,c("fuel","tax","inc","dlic","road")])

 lm.fuel <- lm(fuel ~ tax + inc + dlic + road, data=fuel74)
 coef(lm.fuel)

################# ESIMERKKI 1.3 ############################
 
 attach(fuel74)

 fit <- fitted(lm.fuel)
 res <- resid(lm.fuel)

 plot(fit, res)  # Ei pitäisi näkyä mitään selvää riippuvuutta. Plottaa tämä!
 plot(fuel, res) # Tässä näkyy riippuvuus. Se on matemaattinen fakta. 
                 # Ei ole hyödyllinen kuvio
 plot(fit, fuel); abline(0,1) # Joskus hyödyllinen

 R2 <- sum((fit - mean(fit))^2)/sum((fuel - mean(fuel))^2)

 sigma2.hat <- sum(res^2)/lm.fuel$df.residual
 sigma.hat <- sqrt(sigma2.hat)

 X <- as.matrix(lm.fuel$model)
 X[,1] <- 1
 colnames(X)[1] <- "(Intercept)"
 XX <- t(X)%*%X
 XXinv <- solve(XX)
 coef.se <- sigma.hat*sqrt(diag(XXinv))
 coef.se

 b <- coef(lm.fuel)
 tval <- qt(.975, df=lm.fuel$df.residual)
 conf.int <- rbind(b - tval*coef.se, b + tval*coef.se)
 rownames(conf.int) <- c("2.5%", "97.5%")
 conf.int

 summary(lm.fuel)

################# KUVA 1.3 ##################

 curve(dt(x, df=30), from=-5, to=5, ylab="",xlab="")
 abline(h=0)

 curve(dt(x, df=30, ncp=1), lty=2, from=-5, to=5, ylab="", xlab="",
 add=TRUE)

 curve(dt(x, df=30, ncp=-1), lty=3, from=-5, to=5,
 ylab="", xlab="", add=TRUE)

 q <- qt(.975, df=30)
 lines(c(-q,q),c(dt(-q,df=30,ncp=-1), dt(q,df=30,ncp=1)), type="h")

################## KUVA 1.4 ####################

 curve(1-pt(q,df=30,ncp=x)+pt(-q,df=30,ncp=x), from=-5,to=5,
 ylab="Voimakkuus", xlab="Epäkeskisyysparametri", xaxp=c(-5,5,10),
 yaxp=c(0,1,10))

################## LAPLACE VIRHE S.23 ################

 set.seed(12342214)
 t.sim <- matrix(0,ncol=5,nrow=9999)
 colnames(t.sim) <- names(coef(lm.fuel))
 for (i in 1:9999){
      y <- (log(runif(48)) - log(runif(48)))     # Laplace simulointi
      y.lm <- lm(y~tax+inc+dlic+road, data=fuel74)
      t.sim[i,] <- summary(y.lm)$coef[,"t value"]
                 }

 se.coef <- summary(lm.fuel)$coef[,"Std. Error"]
 coef(lm.fuel)["tax"] +
 sort(t.sim[,"tax"])[c(250,9750)]*se.coef["tax"]

################# PERSENTIILI t BOOTSTRAP S. 24 ##############

 t.boot <- matrix(0,ncol=5,nrow=9999)
 colnames(t.boot) <- names(coef(lm.fuel))
 for (i in 1:9999){
  y <- sample(x=lm.fuel$resid, size=48, replace=TRUE) ## Bootstrap-otos
  y.lm <- lm(y~tax+inc+dlic+road, data=fuel74)
  t.boot[i,] <- summary(y.lm)$coef[,"t value"]
                 }
 coef(lm.fuel)["tax"] + sort(t.boot[,"tax"])[c(250,9750)]*se.coef["tax"]

################ BOOTSTRAP RIVEITTÄIN  S. 24  #################

yX <- as.matrix(fuel74[,c("fuel","tax","inc","dlic","road")])
 boot.coef <- matrix(0,ncol=5,nrow=9999)
 colnames(boot.coef) <- c("(Intercept)","tax","inc","dlic","road")
 for (i in 1:9999){
 s <- sample(1:48, size=48, replace=TRUE)  ## Bootstrap-otosideksit
 y.lm <- lm(yX[s,1]~yX[s,2:5])
 boot.coef[i,] <- coef(y.lm)
                 }

sort(boot.coef[,"tax"])[c(250,9750)]

################# ESIMERKKI 1.5 ###################################

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

yoeng.lm <- lm(formula = PENGL ~ ENGLANTI + RUOTSI + AID1, data = yoeng)

summary(yoeng.lm)

 RSS <- sum(resid(yoeng.lm)^2)
 df <- yoeng.lm$df.resid
 yoeng.lm0 <- lm(PENGL ~ ENGLANTI, data = yoeng)
 RSS0 <- sum(resid(yoeng.lm0)^2)
 df0 = yoeng.lm0$df.resid
 Fobs <- ((RSS0 - RSS)/(df0 - df))/(RSS/df)
 pval <- 1- pf(Fobs, df1=df0-df, df2=df)
 pval

no.eng.lm <- lm(formula = PENGL ~ RUOTSI + AID1, data = yoeng)
summary(no.eng.lm)

################## ESIMERKKI 1.6 ####################################

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

 puut.lm <- lm(formula = vol ~ dbh + ht + d5, data = puut)

 summary(puut.lm)

 uusi <- data.frame(matrix(c(30,50,30,32,25,45), 2,3))
 colnames(uusi) <- colnames(puut)[1:3]
 puut.pred <- predict.lm(puut.lm,newdata=uusi,interval="prediction",
 level=0.9)
 round(puut.pred,2)