Exercise 3.6

Love and money appear to be significant at a 1% level

data(happy, package="faraway")
#str(happy)
lmods <-lm(happy ~ money + sex + love + work, happy)
anova(lmods)
summary(lmods)
## 
## Call:
## lm(formula = happy ~ money + sex + love + work, data = happy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7186 -0.5779 -0.1172  0.6340  2.0651 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.072081   0.852543  -0.085   0.9331    
## money        0.009578   0.005213   1.837   0.0749 .  
## sex         -0.149008   0.418525  -0.356   0.7240    
## love         1.919279   0.295451   6.496 1.97e-07 ***
## work         0.476079   0.199389   2.388   0.0227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.058 on 34 degrees of freedom
## Multiple R-squared:  0.7102, Adjusted R-squared:  0.6761 
## F-statistic: 20.83 on 4 and 34 DF,  p-value: 9.364e-09
confint(lmods, level = .98)
##                      1 %      99 %
## (Intercept) -2.153265677 2.0091037
## money       -0.003147548 0.0223037
## sex         -1.170690142 0.8726746
## love         1.198037838 2.6405203
## work        -0.010659816 0.9628175

** The permutation t-statitics**

nullmod<-lm(happy ~ 1, happy)
anova(nullmod, lmods)
lms <-summary(lmods)

(rss0 <-deviance(nullmod))
## [1] 131.4359
(rss <- deviance(lmods))
## [1] 38.08703
(df0 <- df.residual((nullmod)))
## [1] 38
(df<- df.residual(lmods))
## [1] 34
(fstat<- ((rss0-rss)/(df0-df))/(rss/df))
## [1] 20.83296
1-pf(fstat, df0-df, df)
## [1] 9.364274e-09
nreps<-4000
tstats <- numeric(nreps)
set.seed(123)
for(i in 1:nreps)
{
    lmods <- lm(happy ~ money + sex + love + work, happy)
    tstats[i] <- summary(lmods)$coef[3,3]
}
mean(fstat > lms$fstatistic[1])
## [1] 0

Histogram of test statistics

#d
hist(fstat)

#e

**Boostrap with 90% and 95% confidence levels

 set.seed(123)
 nb<-4000
 coefmat <-matrix(NA,nb,5)
 resids <-residuals(lmods)
 preds<-fitted(lmods)
 for(i in 1:nb){
     booty <- preds + sample(resids, rep=TRUE)
     bmod <- update(lmods, booty ~ .)
     coefmat[i,]<-coef(bmod)
 }
 
colnames(coefmat) <- c("Intercept", colnames(happy[1:4]))
coefmat<-data.frame(coefmat)
apply(coefmat,2,function(x) quantile(x,c(0.10,0.90)))
##      Intercept       happy      money      sex      love
## 10% -1.0816981 0.003456815 -0.6311999 1.556414 0.2407192
## 90%  0.9351802 0.015805139  0.3575996 2.261520 0.7142061
apply(coefmat,2,function(x) quantile(x,c(0.05,0.95)))
##     Intercept      happy      money      sex      love
## 5%  -1.408088 0.00187984 -0.7771200 1.461498 0.1726625
## 95%  1.220105 0.01748788  0.4902119 2.373722 0.7854940