Test for lack of fit

The example can be found in the following link: link

# attach dataset strongx 
data(strongx)

# fit data with simple linear regression model
fit1 <- lm(crossx ~ energy, data=strongx)
summary(fit1)

## 
## Call:
## lm(formula = crossx ~ energy, data = strongx)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.773  -9.319  -2.829   5.571  19.817 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   135.00      10.08    13.4 9.21e-07 ***
## energy        619.71      47.68    13.0 1.16e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.69 on 8 degrees of freedom
## Multiple R-squared:  0.9548, Adjusted R-squared:  0.9491 
## F-statistic: 168.9 on 1 and 8 DF,  p-value: 1.165e-06

par(mfrow = c(2,2))
plot(fit1)

# make a graph of linear regression model
ggplot(data=strongx, aes(y=crossx, x = energy))+
  geom_point()+
  geom_smooth(method="lm", se = FALSE)

Trường hợp này dễ dàng nhận thấy mô hình tuyến tính bậc 1 y=ax+b không phải là mô hình đúng mặc dù giá trị R2 lớn (R2=0.95). Nhìn vào đồ thị Residuals vs fitted thì equal variance không đảm bảo.

data("corrosion")
p <- ggplot(data=corrosion, aes(x=Fe, y = loss))+
  geom_point()+
  geom_smooth(method="lm", se = FALSE)
p

fit2 <- lm(data=corrosion, loss~Fe)
summary(fit2)

## 
## Call:
## lm(formula = loss ~ Fe, data = corrosion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7980 -1.9464  0.2971  0.9924  5.7429 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  129.787      1.403   92.52  < 2e-16 ***
## Fe           -24.020      1.280  -18.77 1.06e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058 on 11 degrees of freedom
## Multiple R-squared:  0.9697, Adjusted R-squared:  0.967 
## F-statistic: 352.3 on 1 and 11 DF,  p-value: 1.055e-09

ga <- lm(loss ~ factor(Fe), data= corrosion)
corrosion <- cbind(corrosion, ga$fit)
p + geom_point(aes(Fe, ga$fit), color="red")

par(mfrow=c(2,2))
plot(fit2)

anova(fit2, ga)

## Analysis of Variance Table
## 
## Model 1: loss ~ Fe
## Model 2: loss ~ factor(Fe)
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1     11 102.850                                
## 2      6  11.782  5    91.069 9.2756 0.008623 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

leveneTest(loss~factor(Fe), data=corrosion)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  6  0.9069 0.5457
##        6

Trường hợp này mặc dù R2=0.97 nhưng mô hình không vượt qua được lack-of-fit test. Nhìn vào đồ thị Residual vs fitted cũng thấy không đảm bảo. Mặc dù Levene test nói ok.

Test for lack of fit

LDA

Jan, 2016