Stat 141  Nov 29
Association among measured variables: Introduction to correlation and regression.

>score12 = read.table(file="D:\\stat141\\th12.dat", header = T)
> plot( 0:25, pch=0:25) # get picture of point options
> plot(th1, th2, pch = 20, main = "Scatterplot: Take Home 1, Take Home 2")
> cor(th1,th2,  use="pairwise" )
[1] 0.5237157

> cor.test(th1,th2,  use="pairwise" )
        Pearson's product-moment correlation
data:  th1 and th2 
t = 5.7997, df = 89, p-value = 9.985e-08
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval: 0.3561880 0.6586236 
sample estimates: cor 0.5237157 

>  cor.test(th1,th2, c = .99, use="pairwise" )
        Pearson's product-moment correlation
data:  th1 and th2 
t = 5.7997, df = 89, p-value = 9.985e-08
alternative hypothesis: true correlation is not equal to 0 
99 percent confidence interval: 0.2975794 0.6942068 
sample estimates:  cor 0.5237157 

> platelet = read.table(file="D:\\stat141\\platelet.dat", header = T) #SW p.550
> #recreate publication Figure in text
> plot(BP, PC, main = "SW Figure 12.14 p.551", xlab = "Blood pressure (mm Hg)", 
           ylab = "Platelet Calcium (nM)", pch =15, bty = "l")

> cor.test(BP,PC, c = .99 )
        Pearson's product-moment correlation
data:  BP and PC 
t = 4.3071, df = 36, p-value = 0.0001219
alternative hypothesis: true correlation is not equal to 0 
99 percent confidence interval: 0.2277740 0.8014416 
sample estimates: cor 0.5831581 

> cor.test(BP,PC )
        Pearson's product-moment correlation
data:  BP and PC 
t = 4.3071, df = 36, p-value = 0.0001219
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval: 0.3238485 0.7609758 
sample estimates: cor 0.5831581
 
> #matches SW ex 12.30 Fishers z-transform, 
   p.563 by hand matches Sw p.561 t-statistic rho = 0

> #correlation measure of linear association SW p566 exs, soybean p.572 ex 12.38
> #Do Soybean at board, illustrate rule of bulge


> ratspeed = read.table(file="D:\\stat141\\ratspeed.dat", header = T) 
                       # fit regression line to amphetamine data ex 12.1
> plot(dose, food, pch = 20)
> tapply( food,dose, mean)
      0     2.5       5 
99.9875 75.5000 54.9500 
> plot(c(0,2.5,5),tapply( food,dose, mean), pch = 15) 
                 #fig 12.10 plot of ave(Y|X) for X = 0, 2.5,5
> ratreg = lm(food ~ dose)  # object for the food on dose straight-line fit
> summary(ratreg)
Call:lm(formula = food ~ dose)
Residuals:
    Min      1Q  Median      3Q     Max 
-21.513  -7.031   1.528   7.448  27.006 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   99.331      3.680   26.99  < 2e-16 ***
dose          -9.008      1.140   -7.90 7.27e-08 ***
---
Residual standard error: 11.4 on 22 degrees of freedom
Multiple R-Squared: 0.7394,     Adjusted R-squared: 0.7275 
F-statistic: 62.41 on 1 and 22 DF,  p-value: 7.265e-08 

> plot(dose, food, pch = 20)
>  abline(ratreg)

> anova(ratreg)
Analysis of Variance Table

Response: food
          Df Sum Sq Mean Sq F value    Pr(>F)    
dose       1 8113.5  8113.5  62.415 7.265e-08 ***
Residuals 22 2859.9   130.0
    
> sqrt(130.0)  [1] 11.40175

> coef(ratreg)
(Intercept)        dose 
   99.33125    -9.00750 
#confidence interval for slope
> coef(ratreg)[2] - qt(.975,22)*1.140;  coef(ratreg)[2] + qt(.975,22)*1.140 
     dose 
-11.37172 
     dose 
-6.643285 

> qqnorm(resid(ratreg))  # not perfect but not awful normality graphic
> plot(fitted(ratreg), residuals(ratreg), pch = 20) #resid vs fit

> predict.lm(ratreg, se.fit=T, interval = "confidence", level=.95)
$fit
        fit      lwr       upr
1  99.33125 91.69979 106.96271
9  76.81250 71.98594  81.63906
17 54.29375 46.66229  61.92521

$se.fit
 [1] 3.679811 3.679811 3.679811 3.679811 3.679811 3.679811
 [7] 3.679811 3.679811 2.327317 2.327317 2.327317 2.327317
[13] 2.327317 2.327317 2.327317 2.327317 3.679811 3.679811
[19] 3.679811 3.679811 3.679811 3.679811 3.679811 3.679811

> predict.lm(ratreg, se.fit=T, interval = "prediction", level=.95)
$fit
        fit      lwr       upr
1  99.33125 74.48502 124.17748
9  76.81250 52.67970 100.94530
17 54.29375 29.44752  79.13998