# multiple regression example
house=read.table("house.txt")
names(house)=c("price","assess","old","months")
pairs(house)
attach(house)

house.lm=lm(price ~ assess + old + months)
summary(house.lm)

plot(assess,price)
lines(assess,fitted(house.lm))    # this doesn't work very well
plot(assess,price)
lines(assess[order(assess)],fitted(house.lm)[order(assess)])  # still weird

plot(fitted(house.lm),resid(house.lm))
plot(assess,resid(house.lm))
plot(months,resid(house.lm))
qqnorm(resid(house.lm))

house.lm2=lm(price ~ assess + I(assess^2) + old + months)
summary(house2.lm)
plot(fitted(house.lm2),resid(house.lm2))
qqnorm(resid(house.lm2))

house.lm3=lm(price ~ assess + I(assess^2) + months)
summary(house.lm3)
plot(fitted(house.lm3),resid(house.lm3))
qqnorm(resid(house.lm3))
predict(house.lm3,data.frame(assess=80,months=6),interval="predict")


# another model building example
oxide=read.table("oxide.txt",skip=18)
names(oxide)=c("day","wind","temp","humid","insol","ox")
attach(oxide)
pairs(oxide)
ox.lm1=lm(ox ~ day + insol + humid + temp + wind)
summary(ox.lm1)
     Coefficients:
     		  Estimate Std. Error t value Pr(>|t|)    
     (Intercept) -12.04010   21.20961  -0.568  0.57553    
     day          -0.02997    0.13995  -0.214  0.83227    
     insol         0.01822    0.05583   0.326  0.74694    
     humid         0.06818    0.13336   0.511  0.61384    
     temp          0.55714    0.15347   3.630  0.00133 ** 
     wind         -0.44749    0.09103  -4.916 5.14e-05 ***
     Residual standard error: 2.977 on 24 degrees of freedom
     Multiple R-Squared: 0.7984,     Adjusted R-squared: 0.7564 
     F-statistic: 19.01 on 5 and 24 DF,  p-value: 1.203e-07 

plot(fitted(ox.lm1),resid(ox.lm1))
pairs(cbind(resid(ox.lm1),oxide[,1:5]))
ox.lm2=lm(ox ~ day + I(day^2) + insol + humid + temp + wind)
summary(ox.lm2)

     (Intercept) -33.410163  18.230624  -1.833  0.07983 .  
     day          -0.730849   0.223740  -3.267  0.00339 ** 
     I(day^2)      0.029055   0.007984   3.639  0.00137 ** 
     insol         0.006761   0.045539   0.148  0.88327    
     humid         0.287954   0.124192   2.319  0.02966 *  
     temp          0.667145   0.128488   5.192 2.90e-05 ***
     wind         -0.406099   0.074940  -5.419 1.66e-05 ***
     Multiple R-Squared: 0.872,      Adjusted R-squared: 0.8387 

pairs(cbind(resid(ox.lm2),fitted(ox.lm2),oxide[,1:5]))
qqnorm(resid(ox.lm2))
ox.lm3=lm(ox ~ day + I(day^2) + humid + temp + wind)
summary(ox.lm3)
     (Intercept) -33.015775  17.664764  -1.869  0.07387 .  
     day          -0.739218   0.212066  -3.486  0.00191 ** 
     I(day^2)      0.029137   0.007801   3.735  0.00103 ** 
     humid         0.284712   0.119740   2.378  0.02573 *  
     temp          0.673103   0.119549   5.630 8.51e-06 ***
     wind         -0.407790   0.072544  -5.621 8.71e-06 ***
     Multiple R-Squared: 0.8719,     Adjusted R-squared: 0.8452 

pairs(cbind(resid(ox.lm3),fitted(ox.lm3),oxide[,1:5]))
qqnorm(resid(ox.lm3))

ox.lm4=lm(ox ~ insol + humid + temp + wind)
summary(ox.lm4)
     (Intercept) -15.49370   13.50647  -1.147  0.26219    
     insol         0.02275    0.05067   0.449  0.65728    
     humid         0.09292    0.06535   1.422  0.16743    
     temp          0.56933    0.13977   4.073  0.00041 ***
     wind         -0.44291    0.08678  -5.104 2.85e-05 ***
     Multiple R-Squared: 0.798,      Adjusted R-squared: 0.7657 

ox.lm5=lm(ox ~ humid + temp + wind)
summary(ox.lm5)
ox.lm6=lm(ox ~ temp + wind)
summary(ox.lm6)
pairs(cbind(resid(ox.lm6),fitted(ox.lm6),temp,wind))
qqnorm(resid(ox.lm6))