STAT 333 Discussion 8 Mar 3/20, 2013Practice Problems1. > x1=c(6,11,14,22,8,15,8,29,14,19,16,7)> x2=c(0.14,0.19,0.23,0.16,0.32,0.27,0.15,0.18,0.22,0.19,0.15,0.24)> x3=c(2.3,4.8,1.7,4.2,3.0,4.7,1.6,2.7,6.1,4.6,1.2,3.8)> x4=c(34,55,32,43,27,58,45,72,39,43,59,30)> y=c(11.6,14.7,9.6,10.3,16.3,18.4,12.0,8.2,17.3,11.1,13.0,12.4)> out1=lm(y~x1+x2+x3+x4)> summary(out1)Call:lm(formula = y ~ x1 + x2 + x3 + x4)Residuals:Min 1Q Median 3Q Max-2.43656 -0.94132 0.01734 0.90497 2.41839Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 3.22524 3.58603 0.899 0.3983x1 -0.35808 0.11470 -3.122 0.0168x2 26.60084 11.83804 2.247 0.0595x3 0.97461 0.39952 2.439 0.0448x4 0.13434 0.05863 2.291 0.0557Residual standard error: 1.898 on 7 degrees of freedomMultiple R-squared: 0.771, Adjusted R-squared: 0.6402F-statistic: 5.893 on 4 and 7 DF, p-value: 0.021252. Use the “bare hand” matrix techqinue to recovery the t table for x1, x2, x3and x4.> X=cbind(rep(1,12),x1,x2,x3,x4);> print(X);x1 x2 x3 x4[1,] 1 6 0.14 2.3 34[2,] 1 11 0.19 4.8 55[3,] 1 14 0.23 1.7 32[4,] 1 22 0.16 4.2 43[5,] 1 8 0.32 3.0 27[6,] 1 15 0.27 4.7 58[7,] 1 8 0.15 1.6 45[8,] 1 29 0.18 2.7 72[9,] 1 14 0.22 6.1 39[10,] 1 19 0.19 4.6 43[11,] 1 16 0.15 1.2 59[12,] 1 7 0.24 3.8 30Now we show:bβ = (X>X)−1X>Y .> betahat=solve( t(X)%*%X ) %*% t(X)%*%y> print(betahat)1[,1]3.2252427x1 -0.3580773x2 26.6008446x3 0.9746084x4 0.1343380What are the estimated standard deviation of each component inbβ? From the result in HW assignment 6,we know that Var(bβ) = σ2(X>X)−1. First, let us check how to estimate σ2. From HW assignment 6, thereis a neat formula for bσ2, which is defined byY>(I−Hx)Yn−p−1.> residual=y %*% (diag(12)-X %*% solve( t(X)%*%X ) %*% t(X) )%*%y/(12-4-1);> print( sqrt(residual) )[,1][1,] 1.898182Then it is straightforward to calculate bσ2(X>X)−1as follows:> cov=solve(t(X)%*%X )*1.898^2;> print( cov )x1 x2 x3 x412.8571197 0.027815595 -32.00344508 -0.1531117181 -0.1323348505x1 0.0278156 0.013153989 0.04042025 -0.0080297383 -0.0043363619x2 -32.0034451 0.040420246 140.11241660 -1.5278636677 0.1816016280x3 -0.1531117 -0.008029738 -1.52786367 0.1595817702 0.0007958579x4 -0.1323349 -0.004336362 0.18160163 0.0007958579 0.0034364320> std=NULL;> for(i in 1:5)+ {+ std[i]=sqrt( cov[i,i] )+ }> print(std)[1] 3.58568259 0.11469084 11.83690908 0.39947687 0.05862109The last question is how to obtain the p-value? This can be completed by> tvalue=betahat/std;> print(tvalue);[,1]0.899478x1 -3.122109x2 2.247280x3 2.439712x4 2.291633The p-values are corresponding to two tail in the t distribution with degree of freedom 7.> pvalue=2*pt(-abs(tvalue),7)> print(pvalue)2[,1]0.39828295x1 0.01679408x2 0.05943615x3 0.04477789x4
View Full Document
Unlocking...