Week 2 RQ 3 get whiteside data > library(MASS) Warning message: package 'MASS' was built under R version 2.2.1 > data(whiteside) > summary(whiteside) Insul Temp Gas Before:26 Min. :-0.800 Min. :1.300 After :30 1st Qu.: 3.050 1st Qu.:3.500 Median : 4.900 Median :3.950 Mean : 4.875 Mean :4.071 3rd Qu.: 7.125 3rd Qu.:4.625 Max. :10.200 Max. :7.200 > attach(whiteside) > # regression that allows separate within-group slopes lesson--have to massage the grouping Insul variable --------------- > InsulN = as.numeric(Insul) > InsulN [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 [39] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 > Insul [1] Before Before Before Before Before Before Before Before Before Before [11] Before Before Before Before Before Before Before Before Before Before [21] Before Before Before Before Before Before After After After After [31] After After After After After After After After After After [41] After After After After After After After After After After [51] After After After After After After Levels: Before After > cnrl = lm(Gas ~ I(InsulN - 1) + Temp + I((InsulN - 1)*Temp)) > summary(cnrl) Call: lm(formula = Gas ~ I(InsulN - 1) + Temp + I((InsulN - 1) * Temp)) Residuals: Min 1Q Median 3Q Max -0.97802 -0.18011 0.03757 0.20930 0.63803 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.85383 0.13596 50.409 < 2e-16 *** I(InsulN - 1) -2.12998 0.18009 -11.827 2.32e-16 *** Temp -0.39324 0.02249 -17.487 < 2e-16 *** I((InsulN - 1) * Temp) 0.11530 0.03211 3.591 0.00073 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.323 on 52 degrees of freedom Multiple R-Squared: 0.9277, Adjusted R-squared: 0.9235 F-statistic: 222.3 on 3 and 52 DF, p-value: < 2.2e-16 > # do within-group regression to check > regB = lm(Gas[Insul =="Before"] ~ Temp[Insul == "Before"]) > summary(regB) Call: lm(formula = Gas[Insul == "Before"] ~ Temp[Insul == "Before"]) Residuals: Min 1Q Median 3Q Max -0.62020 -0.19947 0.06068 0.16770 0.59778 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.85383 0.11842 57.88 <2e-16 *** Temp[Insul == "Before"] -0.39324 0.01959 -20.08 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2813 on 24 degrees of freedom Multiple R-Squared: 0.9438, Adjusted R-squared: 0.9415 F-statistic: 403.1 on 1 and 24 DF, p-value: < 2.2e-16 > regA = lm(Gas[Insul =="After"] ~ Temp[Insul == "After"]) > summary(regA) Call: lm(formula = Gas[Insul == "After"] ~ Temp[Insul == "After"]) Residuals: Min 1Q Median 3Q Max -0.97802 -0.11082 0.02672 0.25294 0.63803 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.72385 0.12974 36.41 < 2e-16 *** Temp[Insul == "After"] -0.27793 0.02518 -11.04 1.05e-11 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3548 on 28 degrees of freedom Multiple R-Squared: 0.8131, Adjusted R-squared: 0.8064 F-statistic: 121.8 on 1 and 28 DF, p-value: 1.046e-11 OK within group regressions check > # significant difference of slopes, do cnrl > vcov(cnrl) (Intercept) I(InsulN - 1) Temp (Intercept) 0.018486202 -0.018486202 -0.0027053168 I(InsulN - 1) -0.018486202 0.032433027 0.0027053168 Temp -0.002705317 0.002705317 0.0005056667 I((InsulN - 1) * Temp) 0.002705317 -0.005050896 -0.0005056667 I((InsulN - 1) * Temp) (Intercept) 0.0027053168 I(InsulN - 1) -0.0050508961 Temp -0.0005056667 I((InsulN - 1) * Temp) 0.0010311886 > 2.13/.1153 #intersection point [1] 18.47355 > Ca = -vcov(cnrl)[2,4]/vcov(cnrl)[4,4] > Ca [1] 4.89813 > DCa = coef(cnrl)[2] + coef(cnrl)[4]*Ca > DCa I(InsulN - 1) -1.565205 > VarDCa = vcov(cnrl)[2,2] + vcov(cnrl)[2,4]*Ca > VarDCa [1] 0.007693079 > # pick-a-point for temp = 4 > diff4 = coef(cnrl)[2] + coef(cnrl)[4]*4 > diff4 I(InsulN - 1) -1.668763 > Vardiff4 = VarDCa + vcov(cnrl)[4,4]*(4 - Ca)^2 > Vardiff4 [1] 0.008524875 > diff4[1] I(InsulN - 1) -1.668763 > diff4[1]/sqrt(Vardiff4) # t-statistic for pick-a-point I(InsulN - 1) -18.07384 > diff4[1] - qt(.975, 52)*sqrt(Vardiff4)[1] # lower endpoint CI I(InsulN - 1) -1.854037 > diff4[1] + qt(.975, 52)*sqrt(Vardiff4)[1] # upper endpoint CI I(InsulN - 1) -1.483488 ---------------- simultaneous region of significance (class handout week 2) do this in R, solve quadratic equation p.318 CNRL paper Also note Week 2 example CNRL data augmented on web page to include full plotting linked in week 2 materials > qf(.95, 2, 52) [1] 3.175141 > cf = coef(cnrl) > vc = vcov(cnrl) > cF = 2*qf(.95, 2, 52) > ?polyroot #polyroot(z) #Arguments: z: the vector of polynomial coefficients in increasing order. > cP = c(cf[2]^2 - cF*vc[2,2], 2*(cf[2]*cf[4] - cF*vc[2,4]), cf[4]^2 - cF*vc[4,4]) > cP I(InsulN - 1) I(InsulN - 1) I((InsulN - 1) * Temp) 4.330847548 -0.427040194 0.006746644 > polyroot(cP) [1] 12.68281+0i 50.61387-0i Simultaneous region: values less than 12.68 and values greater than 50.62 range of Temp measurements approx (0,10) ===================== plotting D(X) and W-H bands # note: need to define n (sample size) in the plot code below; 56 in these data, n-4 is df for the cnrl regression > plot(NA,xlim=c(-3,60),ylim=c(-1,5),xlab="X",ylab="Y",main="Working Hotelling band about D(X)") > abline(h=0) > abline(v=0) > > VarDCa = vcov(cnrl)[2,2] + vcov(cnrl)[2,4]*Ca > lm0 = lm(Gas ~ Temp, subset=(InsulN==1)) > lm1 = lm(Gas ~ Temp, subset=(InsulN==2)) > > abline(lm1$coef[1] - lm0$coef[1], lm1$coef[2] - lm0$coef[2],lwd=3) > temp.seq = seq(-2,25,length=200) > > > dx = (lm1$coef[1] - lm0$coef[1]) + temp.seq*(lm1$coef[2] - lm0$coef[2]) > band = sqrt(2*qf(0.95, 2, n-4) * (VarDCa + vcov(cnrl)[4,4]*(temp.seq - Ca)^2)) > lines(temp.seq, dx - band, lwd=2) > lines(temp.seq, dx + band, lwd=2) > # save your plot window as a pdf before closing with #dev.off() plot at http://statweb.stanford.edu/~rag/stat209/hw5bands.pdf ========================== prior Mathematica copout as a check on cP string get the coeffs for polyroot In[1]:= d[x_] := -2.13 + 0.1153 x 2 vardx[x_] := 0.00769308 + 0.00103119 (x - 4.89813) In[6]:= Expand[vardx[x]] Out[6]= 2 0.032433 - 0.0101018 x + 0.00103119 x In[7]:= Out[4] - 2 3.175 Out[6] Out[7]= 2 4.5369 - 0.491178 x + 0.0132941 x - 2 6.35 (0.032433 - 0.0101018 x + 0.00103119 x ) In[8]:= Expand[%] Out[8]= 2 4.33095 - 0.427032 x + 0.00674604 x > polyroot(c(4.33095, -.427032, .006746)) [1] 12.68321+0i 50.61830-0i --------------------------------------------------------------------------- ------------------------------------------------------------------------