## JOTIMA DEMO 4

#### TEHT 1
airpol<-read.table("http://users.jyu.fi/~junyblom/airpol.dat", header=TRUE)
attach(airpol)
malli<-lm(tmr~perwh+nonpoor+ge65+smean+pmean)
malli0<-lm(tmr~perwh+nonpoor+ge65)
res <- resid(malli)
RSS <- sum(res^2)
res0 <- resid(malli0)
RSS0 <- sum(res0^2)
df <- malli$df
df0 <- malli0$df
OS <-(RSS0-RSS)/(df0-df)
NIM <-RSS/df
F <- OS/NIM
F
p <- 1-pf(F,df1=df0-df,df2=df)
p

# H0 hylätään 

# Toinen tapa:

 anova(malli0, malli, test="F")

#### TEHT 2

### a)
# Sovitteet ja jäännökset
fit <- fitted(malli)
res <- resid(malli)

# Jäänösvarianssi ja -hajonta
sigma2.hat <- sum(res^2)/malli$df.residual
sigma.hat <- sqrt(sigma2.hat)

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

# Luottamusvälit
b <- coef(malli)
tval <- qt(.975, df=malli$df.residual)
conf.int <- rbind(b-tval*coef.se, b+tval*coef.se)
rownames(conf.int) <- c("2.5%", "97.5%")
conf.int

### b)

# t-bootstrap
t.boot <- matrix(0,ncol=6,nrow=9999)
colnames(t.boot) <- names(coef(malli))
for (i in 1:9999){
y <- sample(x=malli$resid, size=50, replace=TRUE) ## Bootstrap-otos
y.lm <- lm(y~perwh+nonpoor+ge65+smean+pmean, data=airpol)
t.boot[i,] <- summary(y.lm)$coef[,"t value"]
}

# Luottamusvälit
se.coef <- summary(malli)$coef[,"Std. Error"]
se.coef
coef(malli)["perwh"] - sort(t.boot[,"perwh"])[c(9750,250)]*se.coef["perwh"]
coef(malli)["nonpoor"] - sort(t.boot[,"nonpoor"])[c(9750,250)]*se.coef["nonpoor"]
coef(malli)["ge65"] - sort(t.boot[,"ge65"])[c(9750,250)]*se.coef["ge65"]
coef(malli)["smean"] - sort(t.boot[,"smean"])[c(9750,250)]*se.coef["smean"]
coef(malli)["pmean"] - sort(t.boot[,"pmean"])[c(9750,250)]*se.coef["pmean"]

# Luottamusvälit lähes samat (hieman leveämmät)

### c)

# Rivi-bootstrap
yX <- as.matrix(airpol[,c("tmr","perwh","nonpoor","ge65","smean","pmean")])
boot.coef <- matrix(0,ncol=6,nrow=9999)
colnames(boot.coef) <- c("(Intercept)","perwh","nonpoor","ge65","smean","pmean")
for (i in 1:9999){
s <- sample(1:50, size=50, replace=TRUE) ## Bootstrap-otosideksit
y.lm <- lm(yX[s,1]~yX[s,2:6])
boot.coef[i,] <- coef(y.lm)
}

# Luottamusvälit
sort(boot.coef[,"perwh"])[c(250,9750)]
sort(boot.coef[,"nonpoor"])[c(250,9750)]
sort(boot.coef[,"ge65"])[c(250,9750)]
sort(boot.coef[,"smean"])[c(250,9750)]
sort(boot.coef[,"pmean"])[c(250,9750)]

# muuttuvat jonkin verran, nämä kannattaisi raportoida


## Teht 3

hooker<-read.table('http://users.jyu.fi/~junyblom/hooker.dat',header=TRUE)
attach(hooker)
Pres<-25.4*Pressure
Tmp<-(5*(Temp-32))/9
logPres<-log(Pres)
malliPres<-lm(logPres~Tmp)
summary(malliPres)
plot(Tmp,logPres)

exp(malliPres$coef[2])
summary(malliPres)

# Pres riippuu kutakuinkin lineaarisesti Tmp:stä
# 1 Celsius-asteen vastaa 3.8% kasvua paineessa


#demo 4 tehtävä 4
osuus<-function(n){
  k<-NULL
  for (i in 1:9999){
    s<-sample(x=1:n,size=n,replace=TRUE)
    z<-1:n %in%s
    k[i]<-mean(z)
    }

  a<-sort(k)
  vali<-c( a[250], a[9750],mean(a) )
 
  }
vali<-osuus(100)
vali
vali<-osuus(1000)
vali
vali<-osuus(5000)
vali

# mitä isompi n sitä kapeampi luottamusväli. 
# keskiarvo pysyy aikalailla samana.