An Introduction Into Linear Statistical Model Building

1 Least Squares In Two Dimensions

###### plot SS 
set.seed(1234)
x    <- seq(0,5,length=30)
y    <- 3 + 5*x + rnorm(length(x),0,1)
x.df <- data.frame(y,x)                   # make a data.frame
res  <- lm(y~x,data=x.df)                 # do linear modeling
summary(res)                              # show results 

## 
## Call:
## lm(formula = y ~ x, data = x.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0306 -0.5751 -0.2109  0.5522  2.7050 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.6800     0.3273   8.188 6.51e-09 ***
## x             5.0094     0.1124  44.562  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9189 on 28 degrees of freedom
## Multiple R-squared:  0.9861, Adjusted R-squared:  0.9856 
## F-statistic:  1986 on 1 and 28 DF,  p-value: < 2.2e-16

So, we know that the p-value < alpha where 2.2e-16 < 0.05, then H0 will be rejected.

a1 <- cov(x,y)/var(x)              # analytical OLS solution 
a0 <- mean(y) - a1*mean(x)
c(a0=a0,a1=a1)                     # show OLS estimators a0,a1

##       a0       a1 
## 2.680009 5.009426

So, the OLS estimators a0 = 2.68 and a1 = 5.01

cor(x.df)                          # shows the Pearson  r.xy

##           y         x
## y 1.0000000 0.9930235
## x 0.9930235 1.0000000

The correlation from x and y is equal to 0.99 or 99%.

sqrt(sum((y - predict(res))^2)/28) # residual standard error

## [1] 0.9188509

cor(x.df$y,predict(res))^2         # R2-value #1: cor(y,y.hat)^2; 

## [1] 0.9860958

ss <- anova(res)
names(ss)                          # what's in ss?

## [1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"

in ss there are Degree of freedoms, Sum of squares, mean of squares, F-value and p-value.

r2 <- ss[,2][1]/sum(ss[,2]) 
r2                                 # R2-value #2: SS.model/ss.total

## [1] 0.9860958

The \(R^2\)-value of the model is 0.98 or 98%. So, the model is good.

F.test <- ss[,3][1]/ss[,3][2]      # Calculate F-value 
F.test

## [1] 1985.773

The F-value is equal to 1985.773

1-pf(F.test,ss[1,1],ss[1,2] )      # Prob(F>F.test)|H0

## [1] 0

The Probability from F-test is 0.