> > #### Fit Ordinary least squares, plot and regressions of dbp, e, |e|, e^2 vs age > > bp.mod1 <- lm(dbp ~ age) > summary(bp.mod1) Call: lm(formula = dbp ~ age) Residuals: Min 1Q Median 3Q Max -16.4786 -5.7877 -0.0784 5.6117 19.7813 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 56.15693 3.99367 14.061 < 2e-16 *** age 0.58003 0.09695 5.983 2.05e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 8.146 on 52 degrees of freedom Multiple R-squared: 0.4077, Adjusted R-squared: 0.3963 F-statistic: 35.79 on 1 and 52 DF, p-value: 2.05e-07 > e <- residuals(bp.mod1) > abs.e <- abs(e) > e2 <- e^2 > > abs_e.mod1 <- lm(abs.e ~ age) # Regression of |e| on age > summary(abs_e.mod1) Call: lm(formula = abs.e ~ age) Residuals: Min 1Q Median 3Q Max -9.7639 -2.7882 -0.1587 3.0757 10.0350 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.54948 2.18692 -0.709 0.48179 age 0.19817 0.05309 3.733 0.00047 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.461 on 52 degrees of freedom Multiple R-squared: 0.2113, Adjusted R-squared: 0.1962 F-statistic: 13.93 on 1 and 52 DF, p-value: 0.0004705 > > e2.mod1 <- lm(e2 ~ age) # Regression of e^2 on age > summary(e2.mod1) Call: lm(formula = e2 ~ age) Residuals: Min 1Q Median 3Q Max -131.59 -46.08 -3.52 31.11 266.55 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -74.2993 36.2691 -2.049 0.045569 * age 3.4921 0.8805 3.966 0.000224 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 73.98 on 52 degrees of freedom Multiple R-squared: 0.2322, Adjusted R-squared: 0.2175 F-statistic: 15.73 on 1 and 52 DF, p-value: 0.0002245 > > par(mfrow=c(2,2)) # Place plots on one page > > plot(age,dbp,main="scatter plot",xlab="age", ylab="dbp") > abline(bp.mod1) > > plot(age,e,main="Residual Plot",xlab="age",ylab="residual") > abline(h=0) > > plot(age,abs.e,main="Absolute Residual Plot",xlab="age",ylab="abs(residual)") > abline(abs_e.mod1) > > plot(age,e2,main="Squared Residual Plot",xlab="age",ylab="residual^2") > abline(e2.mod1) > > s.hat <- predict(abs_e.mod1) # Save the fitted SD's for W > w.hat <- 1/s.hat^2 # Obtain the estimated weights in W > > print(cbind(age,dbp,e,abs.e,s.hat,w.hat)) age dbp e abs.e s.hat w.hat 1 27 73 1.1822391 1.1822391 3.801175 0.069209280 2 21 66 -2.3375761 2.3375761 2.612141 0.146557083 3 22 63 -5.9176069 5.9176069 2.810313 0.126616574 4 24 75 4.9223315 4.9223315 3.206658 0.097251155 5 25 71 0.3423007 0.3423007 3.404830 0.086259931 6 23 70 0.5023623 0.5023623 3.008485 0.110485213 7 20 65 -2.7575453 2.7575453 2.413969 0.171607704 8 20 70 2.2424547 2.2424547 2.413969 0.171607704 9 29 79 6.0221775 6.0221775 4.197519 0.056756367 10 24 72 1.9223315 1.9223315 3.206658 0.097251155 11 25 68 -2.6576993 2.6576993 3.404830 0.086259931 12 28 67 -5.3977917 5.3977917 3.999347 0.062520410 13 26 79 7.7622699 7.7622699 3.603002 0.077031950 14 38 91 12.8019003 12.8019003 5.981070 0.027953887 15 32 76 1.2820851 1.2820851 4.792036 0.043547157 16 33 69 -6.2979457 6.2979457 4.990209 0.040157124 17 31 66 -8.1378841 8.1378841 4.593864 0.047385311 18 34 73 -2.8779765 2.8779765 5.188381 0.037148074 19 37 78 0.3819311 0.3819311 5.782898 0.029902601 20 38 87 8.8019003 8.8019003 5.981070 0.027953887 21 33 76 0.7020543 0.7020543 4.990209 0.040157124 22 35 79 2.5419927 2.5419927 5.386553 0.034464985 23 30 73 -0.5578533 0.5578533 4.395692 0.051754195 24 31 80 5.8621159 5.8621159 4.593864 0.047385311 25 37 68 -9.6180689 9.6180689 5.782898 0.029902601 26 39 75 -3.7781305 3.7781305 6.179242 0.026189640 27 46 89 6.1616539 6.1616539 7.566449 0.017466900 28 49 101 16.4215616 16.4215616 8.160966 0.015014709 29 40 70 -9.3581613 9.3581613 6.377415 0.024587291 30 42 72 -8.5182229 8.5182229 6.773759 0.021794177 31 43 80 -1.0982537 1.0982537 6.971932 0.020572817 32 46 83 0.1616539 0.1616539 7.566449 0.017466900 33 43 75 -6.0982537 6.0982537 6.971932 0.020572817 34 44 71 -10.6782845 10.6782845 7.170104 0.019451321 35 46 80 -2.8383461 2.8383461 7.566449 0.017466900 36 47 96 12.5816232 12.5816232 7.764621 0.016586681 37 45 92 9.7416847 9.7416847 7.368276 0.018419091 38 49 80 -4.5784384 4.5784384 8.160966 0.015014709 39 48 70 -13.9984076 13.9984076 7.962793 0.015771359 40 40 90 10.6418387 10.6418387 6.377415 0.024587291 41 42 85 4.4817771 4.4817771 6.773759 0.021794177 42 55 76 -12.0586232 12.0586232 9.349999 0.011438704 43 54 71 -16.4785924 16.4785924 9.151827 0.011939452 44 57 99 9.7813152 9.7813152 9.746344 0.010527289 45 52 86 -0.3185308 0.3185308 8.755482 0.013044872 46 53 79 -7.8985616 7.8985616 8.953655 0.012473815 47 56 92 3.3613460 3.3613460 9.548172 0.010968811 48 52 85 -1.3185308 1.3185308 8.755482 0.013044872 49 50 71 -14.1584692 14.1584692 8.359138 0.014311232 50 59 90 -0.3787464 0.3787464 10.142689 0.009720617 51 50 91 5.8415308 5.8415308 8.359138 0.014311232 52 52 100 13.6814692 13.6814692 8.755482 0.013044872 53 58 80 -9.7987156 9.7987156 9.944516 0.010111898 54 57 109 19.7813152 19.7813152 9.746344 0.010527289 > > X <- as.matrix(cbind(rep(1,length(dbp)),age)) # X-matrix > Y <- as.matrix(dbp) # Y-vector > W <- as.matrix(diag(w.hat)) # W-hat matrix > > b_w <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% Y # obtain b_w > > sse_w <- t(Y - (X %*% b_w)) %*% W %*% (Y - (X %*% b_w)) # Obtain SSE_w > mse_w <- sse_w/(nrow(X)-ncol(X)) # Obtain MSE_w > > se_bw <- sqrt(diag(mse_w[1,1] * solve(t(X) %*% W %*% X))) # Obtain SE{b_w} > t_bw <- b_w/se_bw # t-statistics > p_bw <- 2*(1-pt(abs(t_bw),nrow(X)-ncol(X))) # P-values > > print(cbind(b_w,se_bw,t_bw,p_bw)) se_bw 55.5657664 2.52091762 22.041881 0.000000e+00 age 0.5963417 0.07923803 7.525953 7.186811e-10 > > ### Fit weighted least squares using lm function and w.hat as weight > > bp.modwls <- lm(dbp ~ age, weight=w.hat) > summary(bp.modwls) Call: lm(formula = dbp ~ age, weights = w.hat) Residuals: Min 1Q Median 3Q Max -2.0230 -0.9939 -0.0327 0.9250 2.2008 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 55.56577 2.52092 22.042 < 2e-16 *** age 0.59634 0.07924 7.526 7.19e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.213 on 52 degrees of freedom Multiple R-squared: 0.5214, Adjusted R-squared: 0.5122 F-statistic: 56.64 on 1 and 52 DF, p-value: 7.187e-10 >