xy = read.table("data030401.txt")

x1 = xy[,1]
x2 = xy[,2]
y = xy[,3];

r1 = lm(y~x1+x2)
summary(r1)

residuals1 = residuals(r1)
par(mfrow = c(2, 2))
plot(x1, residuals1)
plot(x2, residuals1)

 
# there is nonlinearity in the residuals

x11 = x1*x1
x12 = x1*x2
x22 = x2*x2

r2 = lm(y~x1+x2+x11+x12+x22)
# summary(r2)

residuals2 = residuals(r2)
par(mfrow = c(2, 2))
plot(x1, residuals2)
plot(x2, residuals2)


# H0: beta12 = 0, can be accepted

r2a = lm(y~x1+x2+x11+x22)
summary(r2a)

# H0: beta22 = 0, can also be accepted

r2b = lm(y~x1+x2+x11)
summary(r2b)

residuals2b = residuals(r2b)
par(mfrow = c(2, 2))
plot(x1, residuals2b)
plot(x2, residuals2b)

# there is still nonlinearity

x111 = x1^3
x222 = x2^3
x122 = x1*x2*x2
x112 = x1*x1*x2

r3 = lm(y~x1+x2+x11 + x111+ x222 + x122 + x112)
summary(r3)

residuals3 = residuals(r3)
par(mfrow = c(2, 2))
plot(x1, residuals3)
plot(x2, residuals3)


# there is no nonlinearity 

# H0: beta11 = beta222 = beta122 = beta112 = 0

r3a = lm(y~x1+x2 + x111)
summary(r3a)

anova(r3)
anova(r3a)

F = ((0.865-0.836)/4) /(0.836/22) 
F; qf(0.95, 4, 22)

# F = 0.2543860 < F(0.95, 4, 22) = 2.82, accept H0. we accept reduced model, 
# for reduced model

residuals3a = residuals(r3a)
par(mfrow = c(2, 2))
plot(x1, residuals3a)
plot(x2, residuals3a)


# there is no nonlinearity! the final model is 
#  Y = beta0 + beta1*X1 + beta2*X2 + beta111*X^3 + e