# # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- # Ex 13.36 (EPage 543/Anwser EPage 706): # y_hat_fn = function(x1,x2,x3,x4) { 5 + 0.01 * x1 - 0.05 * x2 - 0.13 * x3 - 0.01 * x4 } y_hat = y_hat_fn( 76, 20, 12, 140 ) pnorm( 2.6, mean=y_hat, sd=0.4 ) - pnorm( 1.0, mean=y_hat, sd=0.4 ) # Ex 13.37 (EPage 543/Anwser EPage 706): # y_hat_fn = function(x1,x2) { -0.8 + 0.06 * x1 + 0.9 * x2 } y_hat = y_hat_fn( 50, 3 ) pnorm( 6, mean=y_hat, sd=0.5 ) # Ex 13.38 (EPage 544/Anwser EPage 706): # y_hat_fn = function(x1,x2) { 125.0 + 7.75 * x1 + 0.095 * x2 - 0.0090 * x1 * x2 } y_hat = y_hat_fn( 40, 1100 ) # Ex 13.39 (EPage 544/Anwser EPage 706): # y_hat_fn = function(x1,x2,x3) { 10.0 - 1.2 * x1 + 6.8 * x2 + 15.3 * x3 } y_hat = y_hat_fn( 2, 8, 1 ) y_hat = y_hat_fn( 3, 5, 0 ) # Ex 13.40 (EPage 544/Anwser EPage 706): # y_hat_fn = function(x1,x2,x3,x4) { 1.52 + 0.02 * x1 - 1.4 * x2 + 0.02 * x3 - 0.0006 * x4 } y_hat_fn( 10, 0.5, 50, 100 ) y_hat_fn( 20, 0.5, 10, 30 ) 100 * ( - 0.0006 ) n = 30 SST = 39.2 SSE = 20.0 R2 = 1 - SSE/SST # compute the f statistic: k = 4 numer = R2/k denom = ( 1 - R2 )/( n - (k+1) ) f = numer / denom f_crit = qf( 1-0.05, k, n-(k+1) ) # Ex 13.41 (EPage 544/Anwser EPage 706): # n = 37 R2 = 0.83 k = 6 # compute the f statistic: numer = R2/k denom = ( 1 - R2 )/( n - (k+1) ) f = numer / denom f_crit = qf( 1-0.05, k, n-(k+1) ) # Ex 13.42 (EPage 544/Anwser EPage 706): # beta_2_hat = 3.0 beta_2_se = 0.4321 #t_value = beta_2_hat / beta_2_se t_value = qt( 1-0.05/2, n-(k+1) ) beta_2_ci = beta_2_hat + c(-1,+1) * beta_2_se * t_value y_hat_fn = function(x1,x2) { -200 + 0.21 * x1 + 3.0 * x2 } y_hat = y_hat_fn(1300,7) s_y_hat = 0.353 ci = y_hat + c(-1,+1) * s_y_hat * t_value pi = y_hat + c(-1,+1) * sqrt( s_y_hat^2 + 1.12 ) * t_value # Ex 13.43 (EPage 545/Anwser EPage 706): # DF = read.csv( "../../Data/CH13/ex13-43.txt", header=TRUE, quote="'" ) y_hat_fn = function(x1,x2) { 185.48574 - 45.969466 * x1 - 0.301503 * x2 + 0.088801 * x1 * x2 } y_hat_fn( 2.6, 250 ) 52.0 - y_hat_fn( 2.6, 250 ) y_hat = y_hat_fn( 2.0, 500 ) s_y_hat = 4.69 n = length(DF$y) k = 2 t_value = qt( 1 - 0.05/2, n-(k+1) ) ci = y_hat + c(-1,+1) * s_y_hat * t_value # Ex 13.44 (EPage 545/Anwser EPage 706): # k = 2 n = 25+k+1 # from SPSS output y_hat = 42.253 s_y_hat = 0.350 t_value = qt( 1 - 0.05/2, n-(k+1) ) ci = y_hat + c(-1,+1) * s_y_hat * t_value pi = y_hat + c(-1,+1) * sqrt( s_y_hat^2 + 1.09412^2 ) * t_value t_value = -0.04304 / 0.01773 k = 3 qt( 1 - 0.01/2, n-(k+1) ) # Ex 13.45 (EPage 545/Anwser EPage 706): # n = 25 k = 4 R2 = 0.946 numer = R2/k denom = (1-R2)/(n-(k+1)) f = numer/denom qf( 1-0.05, k, n-(k+1) ) R2_adj = ( (n-1)*R2 - k ) / ( n - (k+1) ) y_hat_fn = function(x1,x2,x3,x4) { 6.121 - 0.082 * x1 + 0.113 * x2 + 0.256 * x3 - 0.219 * x4 } y_hat = y_hat_fn( 16.5, 50, 3, 5 ) y_hat_s = 0.35 t_value = qt( 1 - 0.01/2, n-(k+1) ) pi = y_hat + c(-1,+1) * y_hat_s * t_value # Ex 13.46 (EPage 545/Anwser EPage 706): # n = 12 SSE = 2.09 SST = 12.72 R2 = 1 - SSE/SST k = 2 numer = R2/k denom = (1-R2)/(n-(k+1)) f = numer/denom 1 - pf( f, k, n-(k+1) ) beta_2_hat = 1.250 beta_2_hat_se = 0.312 t_value = beta_2_hat / beta_2_hat_se t_critical = qt( 1-0.05/2, n-(k+1) ) ci = beta_2_hat + c(-1,+1) * beta_2_hat_se * t_critical y_hat_fn = function(x1,x2) { 0.95 + 0.4 * x1 + 1.25 * x2 } y_hat = y_hat_fn( 6, 1 ) y_hat_se = 0.192 sigma2 = SSE/(n-(k+1)) t_value = qt( 1-0.01/2, n-(k+1) ) y_hat + c(-1,+1) * sqrt( y_hat_se^2 + sigma2 ) * t_value # Ex 13.47 (EPage 545/Anwser EPage 706): # n = 30 k = 4 y_hat_fn = function(x1,x2,x3,x4) { 2245 + 28.9 * x1 + 7.64 * x2 + 4.3 * x3 - 37.4 * x4 } y_hat = y_hat_fn( 20, 25, 40, 45 ) y_hat_se = 7.46 t_value = qt( 1-0.05/2, n-(k+1) ) ci = y_hat + c(-1,+1) * y_hat_se * t_value sigma = 31.48 pi = y_hat + c(-1,+1) * sqrt( y_hat_se^2 + sigma^2 ) * t_value # Ex 13.48 (EPage 545/Anwser EPage 706): # n = 15 k = 9 y_hat_fn = function(x1,x2,x3) { 21.967 + 2.8125 * x1 + 1.2750 * x2 + 3.4375 * x3 - 2.208 * x1^2 + 1.867 * x2^2 - 4.208 * x3^2 - 0.975 * x1*x2 - 3.75 * x1*x3 -2.325 * x2*x3 } SSE = 23.379 R2 = 0.938 numer = R2/k denom = (1-R2)/(n-(k+1)) f = numer/denom 1 - pf( f, k, n-(k+1) ) # the p-value for the linear model y_hat = y_hat_fn(0,0,0) y_hat_se = 1.248 t_value = qt( 1-0.05/2, n-(k+1) ) ci = y_hat + c(-1,+1) * y_hat_se * t_value sigma2 = SSE/(n-(k+1)) pi = y_hat + c(-1,+1) * sqrt( y_hat_se^2 + sigma2 ) * t_value l = 3 SSE_full_model = SSE SSE_reduced_model = 203.82 f = ( (SSE_reduced_model - SSE_full_model)/(k-l) ) / ( SSE_full_model/(n-(k+1)) ) f_critical = qf( 1-0.05, k-l, n-(k+1) ) 1 - pf( f, k-l, n-(k+1) ) # the p-value for this test # Ex 13.49 (EPage 545/Anwser EPage 706): # n = 12 k = 2 y_hat_fn = function(x1,x2) { 415.11 - 6.593 * x1 - 4.504 * x2 } y_hat = y_hat_fn(18.9,43) res = 91 - y_hat R2 = 0.768 numer = R2/k denom = (1-R2)/(n-(k+1)) f = numer / denom 1 - pf( f, k, n-(k+1) ) y_hat_se = 8.20 t_value = qt( 1 - 0.05/2, n-(k+1) ) ci = y_hat + c(-1,+1) * y_hat_se * t_value sigma = 24.45 pi = y_hat + c(-1,+1) * sqrt( y_hat_se^2 + sigma^2 ) * t_value y_hat = y_hat_fn(18,45) pi = c(35.94,151.63) pos_negative_se = ( pi - y_hat ) / t_value # this is an interval (-standard_error_y_hat, +standard_error_y_hat) y_hat_se = mean( abs(pos_negative_se) ) # here is the standard error at this input x value t_value = qt( 1 - 0.10/2, n-(k+1) ) # a new t-value pi = y_hat + c(-1,+1) * y_hat_se * t_value # this is a F model utility test beta_1 = 0 vs. beta_1 \ne 0: # DF = read.csv( "../../Data/CH13/ex13-49.txt", header=TRUE, quote="'" ) SST = sum( ( DF$y - mean(DF$y) )^2 ) SSE_full_model = (1-R2) * SST ; k = 2 SSE_reducted_model = (1-0.721) * SST ; l = 1 numer = ( SSE_reducted_model - SSE_full_model )/(k-l) denom = SSE_full_model / ( n-(k+1) ) f = numer / denom (-1.36)^2 # Ex 13.50 (EPage 548/Anwser EPage 706): # t_value = 0.557 / 0.94 k = 5 2*pt( -t_value, n-(k+1) ) # the p-value for this t-value R2 = 0.861 SSE_full_model = (1-R2) * SST ; k = 5 SSE_reducted_model = (1-0.768) * SST ; l = 2 numer = ( SSE_reducted_model - SSE_full_model )/(k-l) denom = SSE_full_model / ( n-(k+1) ) f = numer / denom 1 - pf( f, k-l, n-(k+1) ) # the p-value for this test # Ex 13.51 (EPage 548/Anwser EPage 706): # DF = read.csv( "../../Data/CH13/ex13-51.txt", header=TRUE, quote="'" ) m = lm( y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1*x2), data=DF ) # fit the model suggested summary(m) standardized_residuals = rstandard(m) plot( DF$x1, standardized_residuals ) plot( DF$x2, standardized_residuals ) plot( m$fitted.values, standardized_residuals ) n = length(DF$x1) k = 5 R2 = 0.759 numer = R2/k denom = (1-R2)/(n-(k+1)) f = numer / denom 1-pf( f, k, n-(k+1) ) SST = sum( ( DF$y - mean(DF$y) )^2 ) SSE_full_model = ( 1 - 0.759 ) * SST ; k = 5 SSE_reduced_model = 894.95; l = 2 numer = ( SSE_reduced_model - SSE_full_model ) / ( k-l ) denom = SSE_full_model / ( n-(k+1) ) f = numer / denom 1 - pf( f, k-l, n-(k+1) ) # Ex 13.52 (EPage 549/Anwser EPage 706): # DF = read.csv( "../../Data/CH13/ex13-52.txt", header=TRUE, quote="'" ) n = 20 SST = sum( ( DF$Betacaro - mean(DF$Betacaro) )^2 ) k = 6 + 3 # 6 quadradic terms and 3 linear terms l = 3 SSE_full_model = ( 1 - 0.987 ) * SST SSE_reduced_model = ( 1 - 0.016 ) * SST numer = ( SSE_reduced_model - SSE_full_model ) / ( k-l ) denom = SSE_full_model / ( n-(k+1) ) f = numer / denom 1 - pf( f, k-l, n-(k+1) ) y_hat = 0.66573 y_hat_se = 0.01785 t_value = qt( 1-0.05/2, n-(k+1) ) sigma2 = SSE_full_model / ( n-(k+1) ) pi = y_hat + c(-1,+1) * sqrt( y_hat_se^2 + sigma2 ) * t_value # Ex 13.54 (EPage 549-450/Anwser EPage 706): # SST = 17.2567 n = 31 k = 14 l = 4 SSE_full_model = ( 1 - 0.885 ) * SST SSE_reduced_model = ( 1 - 0.721 ) * SST numer = ( SSE_reduced_model - SSE_full_model ) / ( k-l ) denom = SSE_full_model / ( n-(k+1) ) f = numer / denom 1 - pf( f, k-l, n-(k+1) ) y_hat = 84.5548 y_hat_se = 0.0772 k = 4 t_value = qt( 1 - 0.01/2, n-(k+1) ) ci = y_hat + c(-1,+1) * y_hat_se * t_value