# Question 1


# bee.data <- read.table("D:/Courses/regression/data sets/kaveret.txt",header=F)

num.bees <- bee.data[,2]
hour     <- bee.data[,3]

num.sqrt <- sqrt(num.bees) 

hour.p2 <- hour^2
hour.p3 <- hour^3
hour.p4 <- hour^4
hour.p5 <- hour^5
hour.p6 <- hour^6

lm.3 <- lm(num.sqrt ~ hour + hour.p2)
lm.4 <- lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6)
lm.5 <- lm(num.sqrt ~ factor(hour))

# Alef

anova(lm.3,lm.4)
anova(lm.4)
anova(lm.4)[[2]]
sum(anova(lm.4)[[2]][3:6])


# Bet

x.3 <- model.matrix(lm.3)
x.5 <- model.matrix(lm.5)

dim(x.3) ; x.3[1:10,]
dim(x.5) ; x.5[1:10,]


h.5 <-   x.5 %*% solve(t(x.5)%*%x.5) %*% t(x.5) 

dim(h.5)

lbl <- paste("h-",hour,sep="")
dimnames(h.5) <- list(lbl,lbl)

hour[1:8]
round(h.5[1:8,1:8],4)

table(hour) ; round(1/table(hour),4)

# Gimmel

cbind( model.matrix(lm.3) , h.5 %*%  model.matrix(lm.3))
cbind( model.matrix(lm.4) , h.5 %*%  model.matrix(lm.4))[1:5,]


# Dalet


hour.p7 <- hour^7
hour.p8 <- hour^8
hour.p9 <- hour^9
hour.p10 <- hour^10
hour.p11 <- hour^11
hour.p12 <- hour^12

summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 ))
summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 + hour.p8 ))
summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 + hour.p8 + hour.p9 ))
summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 + hour.p8 + hour.p9 + hour.p10 ))
summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 + hour.p8 + hour.p9 + hour.p10 + hour.p11 ))
summary(lm(num.sqrt ~ hour + hour.p2 + hour.p3 + hour.p4 + hour.p5 + hour.p6 + hour.p7 + hour.p8 + hour.p9 + hour.p10 + hour.p11+ hour.p12))


#   Question 2


# furnace.data <- read.table("E:/Courses/regression/data sets/furnace.txt",header=T)
attach(furnace.data)
ylm <- range(pressure)
xlm <- range(temp)

# Alef 

par(mfcol=c(1,1))
plot(temp[furnace=="A"],pressure[furnace=="A"],pch="A",ylim=ylm,xlim=xlm)
points(temp[furnace=="B"],pressure[furnace=="B"],pch="B")

i.b <- ifelse(furnace=="A",rep(1,19),rep(0,19))
i.a <- ifelse(furnace=="B",rep(1,19),rep(0,19))

temp.a <- ifelse(furnace=="A",temp,rep(0,19))
temp.b <- ifelse(furnace=="B",temp,rep(0,19))

lm.1  <- lm(pressure ~ temp)
lm.1a <- lm(pressure ~ temp,subset=(furnace=="A"))
lm.1b <- lm(pressure ~ temp,subset=(furnace=="B"))

abline(lm.1a,col=2)
abline(lm.1b,col=3)


# Bet 

summary(lm.1)
summary(lm.1a)
summary(lm.1b)


# Gimel + Dalet 

lm.31 <- lm(pressure ~       i.b + temp   + temp.b)
lm.32 <- lm(pressure ~ i.a + i.b + temp.a + temp.b + 0)

summary(lm.31)
summary(lm.32)

summary(lm.1a)
summary(lm.1b)


# Heh + Vav

lm.21 <- lm(pressure ~       temp + temp.b)
lm.22 <- lm(pressure ~ i.b + temp)

summary(lm.22)

anova(lm.21,lm.31)
anova(lm.22,lm.31)

summary(lm.31)
summary(lm(pressure ~ temp * factor(furnace)))


#`	Tesing interaction in lm.32


summary(lm.32)

(est <- rbind(c(0,0,-1,1)) %*% cbind(lm.32$coef))

(se	<-  sqrt(rbind(c(0,0,-1,1)) %*% vcov(lm.32) %*%  cbind(c(0,0,-1,1))))

est / se
2*(1 - pt(est / se,15))