# Here are the class's results from the confidence # interval exercise: > coverage X N_small N_large U_small U_large 1 Jessica Balla 0.95 0.98 0.93 0.95 2 Thais Benoit 0.94 0.98 0.96 0.95 3 Rachel Chan 0.94 0.91 0.95 0.97 4 Jose Godoy 0.98 0.98 0.99 0.97 5 Zoe Griffith 0.96 0.99 0.96 0.93 6 Brendan Harrison 0.96 0.94 0.95 0.93 7 Armin Hojjaty 0.95 0.97 0.94 0.99 8 Zunaira Iqbal 0.91 0.95 0.93 0.95 9 Jonathon Kuntz 0.96 0.95 0.96 0.97 10 Zachary Malone 0.98 0.98 0.96 0.98 11 Jessica Marino 0.93 0.92 0.94 0.94 12 William Meese 0.93 0.93 0.95 0.94 13 Dora Mendez 0.96 0.95 0.97 0.95 14 Danielle Simpson 0.95 0.89 0.94 0.96 15 Jessica Wilson 0.90 0.94 0.90 0.94 > attach(coverage) > stem(N_large) The decimal point is 1 digit(s) to the left of the | 8 | 9 9 | 12344 9 | 555788889 # Each of you, in each condition, performed a series of # 100 independent trials. "Success" occurs if the confidence # interval on one of those trials captures mu. So we have # a series of 100 independent trials with a constant probability # .95 of success (assuming that the confidence intervals are # working right). Those circumstances imply that the distribution # is binomial. # # The most extreme results was 89 successes in the N_large condition. # we can use the binomial distribution to assess how probably that # outcome is: > dbinom(89, 100, .95) [1] 0.007198228 # Not very probable! But keep in mind that across the 15 class members # and 4 conditions we had 60 binomial results. It's not that unexpected # to get a result that should happen only about 1 time in 100 somewhere # across 60 trials. # We can assess how well our results follow a binomial distribution by # doing a binomial Q-Q plot: > temp <- ((1:15)-.5)/15 > Quantiles <- qbinom(temp, 100, .95) > temp [1] 0.03333333 0.10000000 0.16666667 0.23333333 0.30000000 0.36666667 [7] 0.43333333 0.50000000 0.56666667 0.63333333 0.70000000 0.76666667 [13] 0.83333333 0.90000000 0.96666667 > Quantiles [1] 91 92 93 93 94 94 95 95 96 96 96 97 97 98 99 # The plot suggests that the distribution is binomial: > plot(Quantiles, sort(N_large)) # People may have been surprised at the fact that the U3 condition # actually worked pretty well. We can see what the true coverage rate # for that condition is by doing a MUCH larger simulation: > temp <- NULL > for(i in 1:1000000) { + x <- runif(3) + low <- mean(x) -1.96*sqrt(1/12)/sqrt(3) + high <- mean(x) +1.96*sqrt(1/12)/sqrt(3) + temp[i] <- 1 + if(.5 < low) temp[i] <- 0 + if(.5 > high) temp[i] <- 0 + } > mean(temp) [1] 0.952989 # The coverage rate is very close to nominal. So does this mean # that N=3 is sufficiently "large" for the central limit theorem # to work when the population distribution is uniform? # Let's investigate by looking at the sampling distribution of the # mean under those conditions: > temp <- NULL > for(i in 1:1000000) temp[i] <- mean(runif(3)) # The distribution of these means really DOESN'T look normal: > qqnorm(temp); qqline(temp) # The 95% confidence intervals worked pretty well. But will the # coverage rate of, say, 99% intervals also be at the nominal level? > temp <- NULL > for(i in 1:1000000) { + x <- runif(3) + low <- mean(x) -2.576*sqrt(1/12)/sqrt(3) + high <- mean(x) +2.576*sqrt(1/12)/sqrt(3) + temp[i] <- 1 + if(.5 < low) temp[i] <- 0 + if(.5 > high) temp[i] <- 0 + } # No, not so close: > mean(temp) [1] 0.99692 # We continue our investigation of the idea that the # t test is a special form of regression. Note that # the t test in the regression output is identical to # the test reported by the t.test() function: > attach(JackStatlab) > summary(lm(CTRA~CBSEX)) Call: lm(formula = CTRA ~ CBSEX) Residuals: Min 1Q Median 3Q Max -18.9444 -8.0625 0.3715 8.2795 17.0556 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 33.312 1.870 17.817 <2e-16 *** CBSEX -3.368 3.116 -1.081 0.285 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.58 on 48 degrees of freedom Multiple R-squared: 0.02376, Adjusted R-squared: 0.00342 F-statistic: 1.168 on 1 and 48 DF, p-value: 0.2852 > t.test(CTRA~CBSEX, var.equal=TRUE) Two Sample t-test data: CTRA by CBSEX t = 1.0808, df = 48, p-value = 0.2852 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.897595 9.633706 sample estimates: mean in group 0 mean in group 1 33.31250 29.94444 # In fact, most things that we have discussed so far can be # seen as special cases of linear regression. For example, if # we ask R to do a regression with just an intercept and no # slope... > summary(lm(CTRA~1)) Call: lm(formula = CTRA ~ 1) Residuals: Min 1Q Median 3Q Max -21.10 -7.85 0.90 9.40 17.90 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 32.100 1.498 21.42 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.6 on 49 degrees of freedom # ...we get the mean and standard deviation of the variable: > mean(CTRA) [1] 32.1 > sd(CTRA) [1] 10.59505 # Recall that 18 of our children were boys: > sum(CBSEX) [1] 18 # So the first 32 cases are girls: > 50-18 [1] 32 # If we create a variable that has the value -1 for girls and 1 for # boys... > effect <- c(rep(-1,32),rep(1,18)) > effect [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 [26] -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 # ...the regression still does the same t test: > summary(lm(CTRA~effect)) Call: lm(formula = CTRA ~ effect) Residuals: Min 1Q Median 3Q Max -18.9444 -8.0625 0.3715 8.2795 17.0556 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 31.628 1.558 20.299 <2e-16 *** effect -1.684 1.558 -1.081 0.285 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.58 on 48 degrees of freedom Multiple R-squared: 0.02376, Adjusted R-squared: 0.00342 F-statistic: 1.168 on 1 and 48 DF, p-value: 0.2852 # However, this time, the intercept is the mean of the boys' # and girls' means, and the slope is the difference from that # "mean of means" associated with being in one group: > mean(tapply(CTRA, CBSEX, mean)) [1] 31.62847 # That coding approach is called "effects coding." If the # sample sizes are the same for the two groups, the intercept # is the grand (i.e., unconditional) mean, and the slope # represents the "effect" of being in one group or the other. # The fact is, ANY coding system that uses a single variable # to identify the two groups will accomplish the same test: > nonsense <- c(rep(pi, 32),rep(62485, 18)) > nonsense [1] 3.141593 3.141593 3.141593 3.141593 3.141593 [6] 3.141593 3.141593 3.141593 3.141593 3.141593 [11] 3.141593 3.141593 3.141593 3.141593 3.141593 [16] 3.141593 3.141593 3.141593 3.141593 3.141593 [21] 3.141593 3.141593 3.141593 3.141593 3.141593 [26] 3.141593 3.141593 3.141593 3.141593 3.141593 [31] 3.141593 3.141593 62485.000000 62485.000000 62485.000000 [36] 62485.000000 62485.000000 62485.000000 62485.000000 62485.000000 [41] 62485.000000 62485.000000 62485.000000 62485.000000 62485.000000 [46] 62485.000000 62485.000000 62485.000000 62485.000000 62485.000000 > summary(lm(CTRA~nonsense)) Call: lm(formula = CTRA ~ nonsense) Residuals: Min 1Q Median 3Q Max -18.9444 -8.0625 0.3715 8.2795 17.0556 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.331e+01 1.870e+00 17.816 <2e-16 *** nonsense -5.390e-05 4.987e-05 -1.081 0.285 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.58 on 48 degrees of freedom Multiple R-squared: 0.02376, Adjusted R-squared: 0.00342 F-statistic: 1.168 on 1 and 48 DF, p-value: 0.2852 > detach(JackStatlab) # We began our investigation of multiple regression. First, we # looked at the one circumstance under which simple regression # slopes will be identical to multiple regression slopes. # Here's the multiple regression of y on x1 and x2: > attach(uncor) > lm(y~x1+x2) Call: lm(formula = y ~ x1 + x2) Coefficients: (Intercept) x1 x2 10 10 -5 # And here, we see that the simple regression slopes for # x1 and x2 are the same as the slopes in the multiple regression: > lm(y~x1) Call: lm(formula = y ~ x1) Coefficients: (Intercept) x1 -17.5 10.0 > lm(y~x2) Call: lm(formula = y ~ x2) Coefficients: (Intercept) x2 60 -5 # That will happen ONLY if x1 and x2 are perfectly uncorrelated: > cor(x1,x2) [1] 0 # Here's the data set. For the moment, this result is an interesting # bit of trivia, but the idea will turn out to be profoundly useful # later on. > uncor y x1 x2 1 15 1 1 2 -30 1 10 3 30 3 2 4 -5 3 9 5 45 5 3 6 20 5 8 7 60 7 4 8 45 7 7 9 75 9 5 10 70 9 6 > > > # Here's a sample data set for multiple regression. We will # attempt to understand differences in the conditional mean # of a Neuroticism scale using Agreeableness and Depression # as predictors. Here's the data set: > attach(BigFive) > BigFive Neuroticism Agreeableness Depression 1 66 79 15 2 11 79 6 3 7 87 15 4 71 17 17 5 87 79 33 6 60 92 11 7 94 27 32 8 1 79 4 9 32 63 30 10 66 44 8 11 43 87 11 12 55 57 40 13 60 63 7 14 49 50 14 15 27 57 4 16 7 79 1 17 3 87 6 18 96 6 51 19 60 74 12 20 14 57 12 21 9 63 4 22 55 63 6 23 43 63 8 24 32 74 3 25 60 65 14 26 90 17 38 27 60 63 14 28 11 87 4 29 94 69 45 30 76 83 25 31 66 50 3 32 18 90 8 33 37 44 27 34 97 44 13 35 60 63 28 36 18 96 24 37 60 44 11 38 37 69 12 39 55 44 12 40 71 50 31 41 11 27 13 # Getting multiple regression results in R is as simple # as adding predictors to the list: > lm(Neuroticism ~ Depression + Agreeableness) Call: lm(formula = Neuroticism ~ Depression + Agreeableness) Coefficients: (Intercept) Depression Agreeableness 55.1270 1.0918 -0.4049 # But what's really going on? Here, we regress Neuroticism # on Depression (ignoring Agreeableness): > NonD <- lm(Neuroticism~Depression) > NonD Call: lm(formula = Neuroticism ~ Depression) Coefficients: (Intercept) Depression 25.610 1.368 # The residuals from that regression represent variation # in Neuroticism that we CAN'T understand using Depression # as a predictor: > Yres <- NonD$res # Now we regress Agreeableness on Depression: > AonD <- lm(Agreeableness~Depression) > AonD Call: lm(formula = Agreeableness ~ Depression) Coefficients: (Intercept) Depression 72.8911 -0.6809 # The residuals from that regression represent variation in # Agreeableness that doesn't overlap with variation in Depression: > Xres <- AonD$res # So if we look at the relationship between those two sets of # residuals, we are looking at how unique variation in Agreeableness # can help us understand Neuroticism variation that Depression can't # help us understand. (Wow, that was a mouthful!) > lm(Yres~Xres) Call: lm(formula = Yres ~ Xres) Coefficients: (Intercept) Xres -1.055e-16 -4.049e-01 > # Note that the slope from that regression of residuals on residuals # is identical to the slope for Agreeableness in the multiple regression: > lm(Neuroticism ~ Depression + Agreeableness) Call: lm(formula = Neuroticism ~ Depression + Agreeableness) Coefficients: (Intercept) Depression Agreeableness 55.1270 1.0918 -0.4049