Homework Exercises and Solutions Guide

Initializing RStudio

The code block below uses the library function to ensure that the Mosaic package is loaded and will import the data frames we have used in all modules: Data3350, Dolphin, Med, Perc, PercB. Please make sure all data sets are downloaded and placed in your working directory.

library(mosaic)
library(readxl)
Data3350 = read_excel("Data3350.xlsx")
Med = read_excel("Med.xlsx")
Perc = read_excel("Perc.xlsx")
PerB = read_excel("PercB.xlsx")

Homework Set 1: Modules 1 - 3

Item 1.1

Using the built-in R data frame ChickWeight which examines weight vs age of chicks on different diets, conduct an analysis of the variable weight. Include all relevant plots, and include titles and axis labels for your histogram and density plot.

favstats(~weight , data = ChickWeight)

The descriptives indicate possible right skew as distance between Q3 and Max (210) is much greater than distance between Q1 and min (28). Also, by Robb’s Rule of Thumb, the mean is more than 18 units greater than the median which is far more than one-tenth of a standard deviation: \[\frac{s}{10}=\frac{71.07}{10}= 7.11\] Therefore we expect skew and outliers to the right.

We can begin our check for outliers by calculating 2 standard deviations away from the mean. Admittedly, this data set has 578 observations, so we definitely expect several data points to exceed this distance from the mean. (See Module S3 material to understand why I make this claim.)

histogram(~weight , data = ChickWeight)

densityplot(~weight , data = ChickWeight)

Both histogram and density plot have some evidence of a mode (tallest bar or bars) with a tail in both directions. The left tail is severely truncated. Hence, the data appear to have been drawn from an approximately normal distribution, with strong skew to the right.

boxplot(ChickWeight$weight, horizontal = TRUE)

There are 7 outliers to the right, and none to the left, confirming the strong skew to the right which we already suspected.

stem(ChickWeight$weight, scale = 2)

  The decimal point is 1 digit(s) to the right of the |

   3 | 599999999
   4 | 00000111111111111111111112222222222222223333456678888888899999999999
   5 | 000000111111112222333334445555566667778888899999
   6 | 001111111222222223333344444555556666777778888889
   7 | 0011111122222233333444444446667778889999
   8 | 0011222334444445555556677778899999
   9 | 0001223333566666788888889
  10 | 00001111222333333345666677788899
  11 | 01122223445555667789
  12 | 0000222333334444555566778889
  13 | 0113444555566788889
  14 | 1112344445555666667778889
  15 | 0011234444555666777777789
  16 | 0000223333444446678899
  17 | 0000134445555789
  18 | 1224444455567778
  19 | 2225677778889999
  20 | 01234445555579
  21 | 00245578
  22 | 00123577
  23 | 01123344556788
  24 | 08
  25 | 001699
  26 | 12344569
  27 | 259
  28 | 0178
  29 | 0145
  30 | 35579
  31 | 8
  32 | 127
  33 | 12
  34 | 1
  35 | 
  36 | 1
  37 | 3

The stem plot is optional as the sample size makes it unweildy. A scale of 1 is not detailed enough, and a scale of 5 is too granular. A scale of 2 (splitting stems two ways) shows the distribution rather well. We can also pick off the highest 7 observations which were identified as outliers, the values \(x = 32.2, 32.7, 33.1, 33.2,34.1, 36.1, 37.3\).

Summary of EDA. Approximately normal but skewed right distribution with outliers to the right.

Item 1.3

Use the built-in R data frame cars which compares the numeric variables speed measured in miles per hours (mph) and stopping distance (variable is dist, measured in feet). The data is from the 1920’s, by the way. Create an xyplot for the two variables, and include axis labels on your plot.

xyplot(speed ~ dist, data = cars, 
       xlab = "Stopping Distance (ft)",
       ylab = "Speed (mph)")

Item 1.5

Using the built-in R data frame Dimes which compares the mass of a dime to the year in which is was minted, conduct an analysis of the variable mass. Include all relevant plots, and include titles and axis labels for your histogram and density plot.

favstats(~mass, data = Dimes)

The mean and the median differ by approximately a tenth of a standard deviation, we expect skew and outliers to the right since the mean is to the right of the median. However, the tails of the distribution appear to be approximately symmetric. Notice that none of the observations are more than 2 standard deviations from the mean.

histogram(~mass, data = Dimes,
          fit = "normal")

histogram(~mass, data = Dimes,
          fit = "normal",
          width = .01)

The width of the bins is hard to choose, here. The standard histogram with bin width of 0.02 has too few bars to see the distribution, while a bin width of 0.01 is too granular. Let’s check the other graphics before deciding on any analysis statements.

densityplot(~mass, data = Dimes)

The plateau-shaped density plot shows the data appear more uniformly distributed than normally distributed.

boxplot(Dimes$mass, horizontal = TRUE)

We see the approximate symmetry with no outliers.

stem(Dimes$mass)

  The decimal point is 2 digit(s) to the left of the |

  221 | 4
  222 | 
  223 | 03445568
  224 | 79
  225 | 244569
  226 | 888
  227 | 14457
  228 | 27
  229 | 288

The best evidence yet of a bell-shaped distribution. Still, one is surely uncomfortable with the normality assumption given this data. For robust procedures like the \(t\)-test and \(t\)-intervals, we would proceed using the normality since there are no outliers. Outliers cause problems in statistics procedures with small samples sizes because they are influential data points. Here, we have minor issues deciding whether the shape is more uniform than normal, but that won’t affect the accuracy of \(t\) procedures.

Summary of EDA. Approximately normal and approximately symmetric. Minor issues with analysis, but not enough to be troublesome in practice.

Item 2.1

Using the Mosaic function rflip, flip 16 coins and count the number of Heads. Repeat 10,000 times using the Mosaic function do. Estimate the probability that, if Dr. Bristol tasted 16 cups of tea each of which was randomly chosen to be “tea first” or “milk first,” that she would get at least 14 correct using a histogram with the type parameter set to “count.”

coins = do(10000) * rflip(16)
tally(~ heads, data = coins)

heads
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
   4   16   87  296  629 1215 1687 1973 1806 1204  707  270   88   13    4 
  16 
   1 

histogram(~heads, data = coins,
          width = 1,
          type = "count")

We see that 19 times there were exactly 14 Heads, 4 times there were exactly 15 Heads and 1 time there were exactly 16. Thus, we estimate the probability: \[P(k\geq 14)=\frac{24}{10000} = 0.0024\] Notice that, as the sample size grows, Dr. Bristol’s accuracy rate can be much lower and still provide strong evidence to reject the null as unlikely to be true.

Item 2.3

For 32 coins flips (and 10,000 randomized draws), estimate the probability that Dr. Bristol would get at least 24 correct if she were guessing at random. Use a histogram with the type parameter set to “count.” Update pFun to find the theoretical probabilities and compare the two results.

coins = do(10000) * rflip(32)
tally(~ heads, data = coins)

heads
   5    6    7    8    9   10   11   12   13   14   15   16   17   18   19 
   1    2   12   18   52  122  310  502  786 1106 1344 1422 1332 1091  810 
  20   21   22   23   24   25   26   27 
 504  329  162   62   23    5    4    1 

histogram(~heads, data = coins,
          width = 1,
          type = "count")

Notice that we have \[P(k\geq 24)=\frac{30+7+2+1}{10000}=\frac{40}{10000}= 0.004\]

For the pFun function, we update it as follows:

pFun <- function(t) {
  choose(32,t)*(.5)^t *(.5)^(32-t)
}

Sum the probabilities from \(k=24\) to \(k=32\).

sum(pFun(24:32))

[1] 0.003500183

Note that the two values are reasonably close to the same value. We could certainly create more than 10,000 samples to improve the accuracy of the randomization effort at estimating the theoretical probability.

Item 3.1

Find the percentile ranking for an IQ score of 120 using the xpnorm function. IQ’s have the \(N(100,15)\) distribution.

xpnorm(120,100,15)

If X ~ N(100, 15), then 

    P(X <= 120) = P(Z <= 1.333) = 0.9088
    P(X >  120) = P(Z >  1.333) = 0.09121

[1] 0.9087888

An IQ score of 120 corresponds to the 90th percentile - note it’s not the 91st percentile, we don’t round, we truncate.

Item 3.4

Find the approximate SAT-Math score that corresponds to the 60th percentile SAT-Math score. Use the xqnorm function. SAT components have the \(N(500,100)\) distribution.

xqnorm(.6,500,100)

If X ~ N(500, 100), then 

    P(X <= 525.3347) = 0.6
    P(X >  525.3347) = 0.4

[1] 525.3347

The 60th percentile on an SAT component is approximately 525.

Homework Set 2: Modules 4 and 5

Item 4.3

Using the AccDate variable from the Data3350 data frame, test the hypothesis that the Yes/No responses to the dating question depend upon whether or not the person is a member of social Greek fraternity or sorority. Test at the \(\alpha = 0.1\) level, and include a Mosaic plot with a description about it’s relationship to your \(p\)-value and conclusions.

Solution. The data pass verification, and we fail to reject the null. There is no evidence that Yes/No responses to the Dating Question depend upon membership in a social Greek fraternitry or sorority. Discussion will vary.

prop.test(AccDate ~ Greek, data = Data3350, alternative = "greater")

    2-sample test for equality of proportions with continuity correction

data:  tally(AccDate ~ Greek)
X-squared = 0.00045324, df = 1, p-value = 0.4915
alternative hypothesis: greater
95 percent confidence interval:
 -0.1333805  1.0000000
sample estimates:
   prop 1    prop 2 
0.4414414 0.4259259 

mosaicplot(AccDate ~ Greek, data = Data3350)

Item 4.4

A strong sense of Coping Humor indicates a person who uses humor to relieve stress and deal with the struggles of life. Test at the .05 level whether more than 10% of North Georgia students exhibit strong Coping Humor. A recent study used a criteria of scoring 30 or higher on the Coping Humor Scale to evaluate this criteria, and found that 21 of 175 North Georgia students did so.

We are testing the hypothesis \[H_0 : prop = 0.1 \\ H_a : prop > 0.1\] using a 1-proportion \(z\)-test which RStudio will evaluate using \(\chi^2\). The data pass verification for \(z\)-procedures because the estimated number of successes is \[np_0 = (175)(0.1) = 17.5\] and thus, the estimated number of failures is much, much larger than 10.

prop.test( x = 21 , n = 175 , p = 0.1 , 
           alternative = "g",
           correct = F)

    1-sample proportions test without continuity correction

data:  21 out of 175
X-squared = 0.77778, df = 1, p-value = 0.1889
alternative hypothesis: true p is greater than 0.1
95 percent confidence interval:
 0.08527351 1.00000000
sample estimates:
   p 
0.12 

I showed the abbreviation format for the alternative and correct parameters. We fail to reject the null because \(p = 0.19 > 0.05 = \alpha\). There is no evidence that more than 10% of UNG students exhibit strong coping humor.

Item 4.5

A recent study asked UNG students whether they frequently texted at work about things unrelated to work. For younger students, 25 of 103 said they frequently did so while 15 of 46 students who 21 or older reported doing so. Test for an age-difference for Texting Frequently at Work at the .05 level.

If we use the subscripts Y and D for “younger” and “older” students respectively, we are testing the hypothesis \[H_0 : p_Y = p_D \\ H_a : p_Y > p_D\] using a 2-proportion \(z\)-test (calculator) which RStudio will evaluate using \(\chi^2\). The data pass verification for \(z\)-procedures as there more than 10 successes and 10 failures in both samples. The data pass verification for \(\chi^2\) procedures because they are less strict than for \(z\).

As a fun option, we can create a mosaic plot for the problem since it can be represented as a two-way table. You will not need to do this in homework, tests or projects, but it is not difficult. I just like mosaic plots, so I wanted to learn how to create one from summary data alone. We first create a table called tab, then ask for mosaic plot.

tab = matrix(c(25,15,78,31),ncol=2,byrow=TRUE)
colnames(tab) = c("Younger", "Older")
rownames(tab) = c("Yes","No")
tab

    Younger Older
Yes      25    15
No       78    31

The mosaic plot will display the relative proportions as areas.

mosaicplot(tab)

We need to specify a one-tailed test and provide the summary statistics.

prop.test( x = c(25, 15), n = c(103,46), 
           alternative = "greater" ,
           correct = FALSE  )

    2-sample test for equality of proportions without continuity
    correction

data:  c out of c25 out of 10315 out of 46
X-squared = 1.1254, df = 1, p-value = 0.8556
alternative hypothesis: greater
95 percent confidence interval:
 -0.2166097  1.0000000
sample estimates:
   prop 1    prop 2 
0.2427184 0.3260870 

With \(p = .86 > 0.05 = \alpha\), in fact, much greater, we fail to reject the null. There is no evidence younger students text more frequently at work about things not related to work.

Item 5.2

Use the toxic relationship beliefs variable TxRel from the Data3350 data frame to test for a significant difference between females and males (Sex variable). Use the 0.05 level of significance.

This calls for an independent samples t-test. \[H_0 : \mu_F = \mu_M\\H_a : \mu_F \neq \mu_M\] For some reason, the variable TxRel is not being recognized by the Mosaic function histogram. I use the R standard hist, instead, with the traditional R structure of DataFrame$variable. Notice that other functions are working fine.

favstats(TxRel ~ Sex, data = Data3350)

boxplot(Data3350$TxRel, horizontal = TRUE)

hist(Data3350$TxRel)

histogram(TxRel ~ Sex, data = Data3350)

Error in r[i1] - r[-length(r):-(length(r) - lag + 1L)] : 
  non-numeric argument to binary operator

Again, I can’t get this work. I don’t know why. Mosaic’s densityplot and favstats functions will not work either for the variable by itself. Still, we have evidence of a bell-shaped distribution with no outliers and a huge overall sample size, all of which means that the data are quite suitable for \(t\)-procedures. Since the sample sizes are not sharply unequal (ratio of 96:68 is less than 2:1), these data easily pass verification.

Because \[p = 0.01616 < 0.05 = \alpha\] we reject the null. Evidence suggests there is a difference in naive relationship beliefs based on biological sex. A quick glance at the boxplot using the grouping variable Sex clarifies the situation, as will a quick glance at the means summary in the last two lines of the output.

Item 5.3

Use the Anxiety variable Anx from the Data3350 data frame to test for a significant difference between females and males (Sex variable). Use the 0.05 level of significance.

Another independent samples t-test. \[H_0 : \mu_F = \mu_M\\H_a : \mu_F \neq \mu_M\] Even though the sample size is much greater than 40 and easily passes verification, there’s no reason not to check the distribution and descriptives.

We have evidence of an approximately bell-shaped distribution that is skewed right with an outlier to the right. However, outliers only cause inaccurate \(p\)-values for \(t\)-tests when the sample size is less than 40. Here, we have \(n=144\). There is also no issue with the homogeneity of variances assumption as our sample sizes are not sharply unequal. We have 84 females and 60 males which is not greater than a ratio of 2:1.

We reject the null at the 0.05 level and note that would also have rejected the null at the 0.01 level. We have strong evidence that anxiety levels differ based upon biological sex.

Item 5.5

Use the Neuroticism variable Neuro from the Data3350 data frame to test for a significance difference in levels of Neuroticism based upon Primary Humor Style (PHS variable). Conduct a post hoc Tukey HSD test if needed.

Because the grouping variable has four levels, we will use ANOVA.

\[H_0 : \mu_{AF} = \mu_{SE} = \mu_{AG} = \mu_{SE} \\H_a : \text{At least one significantly different}\]

Let’s take a look at the boxplots and descriptive using the grouping variable to get breakdowns for the individual samples. Again, this is not really necessary due to sample size and the robustness of the \(F\) statistic.

The main reason for running the favstats function is to check the relative sample sizes to assess the homegeneity assumption. Notice the largest to smallest group size ratio is 43 to 31 which is far less than 2:1. The sample size is more than adequate, and the group sizes are not sharply unequal. OK to proceed with ANOVA procedure.

modA = lm(Neuro ~ PHS, data = Data3350)
anova(modA)

Analysis of Variance Table

Response: Neuro
           Df  Sum Sq Mean Sq F value    Pr(>F)    
PHS         3  5145.1 1715.03   10.73 2.158e-06 ***
Residuals 140 22375.9  159.83                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

With \[p = 0.00000216 < 0.05 = \alpha\] we reject the null and have evidence for a difference in Neuroticism based on primary humor style. When we reject the null with an ANOVA procedure, we must follow up with a post hoc test. In this course, we will always use Tukey’s Honestly Significant Difference (HSD).

TukeyHSD(modA)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = x)

$PHS
            diff         lwr       upr     p adj
AG-AF   4.978245  -2.7668457 12.723335 0.3427512
SD-AF  12.000872   3.9970309 20.004713 0.0008530
SE-AF  -4.282502 -12.5045100  3.939505 0.5301145
SD-AG   7.022627  -0.3485287 14.393783 0.0679484
SE-AG  -9.260747 -16.8682354 -1.653259 0.0101553
SE-SD -16.283374 -24.1541382 -8.412610 0.0000018

mplot(TukeyHSD(modA))

While we can use the HSD confidence intervals and corresponding \(p\)-values to find the significant pairwise differences, the mplot visual output makes life much easier. Only intervals that do not include zero are significantly different, so we see significant differences in

• Self-Defeating vs. Affiliative
• Self-Enhancing vs. Aggressive
• Self-Enhancing vs. Self-Defeating

Item 5.6

Use the Self-Esteem variable SE from the Data3350 data frame to test for a significance difference in levels of Self-Esteem based upon Primary Humor Style (PHS variable). Conduct a post hoc Tukey HSD test if needed.

We will test the hypothesis

\[H_0 : \mu_{AF} = \mu_{SE} = \mu_{AG} = \mu_{SE} \\H_a : \text{At least one significantly different}\]

using ANOVA.

To check the assumptions, consider:

favstats(SE ~ PHS , data = Data3350)

The overall sample size is far greater than 20, and the largest to smallest group size ratio is less than 2:1, so these data are appropriate for ANOVA procedures. For fun, let’s look at the boxplots by group.

boxplot(SE ~ PHS , data = Data3350)

Now we conduct the test by first running the lm (linear models) procedure and next by calling anova. I am creating the variable modS (for self-esteem model) to store the model to used later.

modS = lm(SE ~ PHS , data = Data3350)
anova(modS)

Analysis of Variance Table

Response: SE
           Df Sum Sq Mean Sq F value  Pr(>F)    
PHS         3   3644 1214.68  10.884 1.8e-06 ***
Residuals 140  15625  111.61                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Since we reject the null and find significant evidence of differences, we must conduct the post hoc Tukey’s HSD procedure.

TukeyHSD(modS)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = x)

$PHS
           diff        lwr        upr     p adj
AG-AF -3.972243 -10.444327  2.4998406 0.3843188
SD-AF -9.646033 -16.334338 -2.9577279 0.0014582
SE-AF  4.182796  -2.687818 11.0534090 0.3916473
SD-AG -5.673790 -11.833400  0.4858203 0.0827024
SE-AG  8.155039   1.797940 14.5121372 0.0059398
SE-SD 13.828829   7.251728 20.4059300 0.0000012

mplot(TukeyHSD(modS))

The two significant pairwise differences are between Self-Enhancing vs. Aggressive and between Self-Enhancing vs. Self-Defeating humor styles.

Item 5.8

Use the subset command as shown above for the following hypothesis test. Males are more likely to use Aggressive Humor, and younger folks are more likely to use aggressive humor than older folks. Find the overall aggressive humor average using the (HSAG variable) in the Data3350 data frame. Subset a group of young men who are less than 20 years old, and test the hypothesis that this subpopulation has a higher group mean than the overall population mean.

First, we must determine \(\mu_0\), the overall mean.

favstats(~HSAG , data = Data3350)

Next, we need a subset of young males.

male = subset(Data3350, Sex == "M", c(Age,HSAG))
yMale = subset(male, Age < 20, HSAG)
yMale

We will test the hypothesis \[H_0 : \mu = 28.85 \\ H_a : \mu > 28.85 \] using a one-tailed \(t\)-test. Checking the assumptions, we need to know our sample size.

favstats(~HSAG , data = yMale)

With a sample size of 12, we have some work to do on the normality assumption: both verify using a boxplot there are no outliers, and use a histogram and density plot to verify the samples appears drawn from an approximately normal distribution.

boxplot(yMale$HSAG, horizontal = T)

We find no outliers, so we check a histogram and density plot. Don’t worry about the extra code - just a demo on how to plot multiple graphics in one output window. It’s totally optional. If you’re interested, though, the output window has coordinates like the unit square, and we define the bottom-left point and top-right point for each graphic inside the print command.

p1 = histogram(~HSAG , data = yMale, width = 3)
p2 = densityplot(~HSAG , data = yMale)
print(p1, c(0,0.5,1,1), more = TRUE)
print(p2, c(0,0,1,0.5))

I had to fiddle a bit with the width parameter on the histogram, but the overall shape confirms the sample appear to have been drawn from an approximately normal distribution. The density plot also shows an approximate bell shape.

t.test(~HSAG , data = yMale , mu = 28.85 , alternative = "g")

    One Sample t-test

data:  HSAG
t = 1.9429, df = 11, p-value = 0.03903
alternative hypothesis: true mean is greater than 28.85
95 percent confidence interval:
 29.1072     Inf
sample estimates:
mean of x 
    32.25 

We reject the null because \(p=0.039 < 0.05 = \alpha\), and we thus have evidence that younger males use more aggressive humor than the overall population at UNG.

Homework Set 3: Modules 6 and 7

Item 6.1

How many resamplings should we use when bootstrapping? Try re-running the code blocks from the Sleep example with 50, 100, 500, and 1000 resamplings. How does the accuracy compare to the theoretical confidence interval as number of resamplings increases? Explain why about 500 resamplings is usually good enough.

Notice how the endpoints fluctuate dramatically as your retry this code block. In five iterations, I saw left endpoints of 6.06, 6.07, 6.09 , 6.08 and 6.02.

bootstrap = do(50) * mean(resample(Data3350$Sleep))
qdata(~mean, p=c(0.025, 0.975), data=bootstrap)

    2.5%    97.5% 
6.018258 6.598182 

Still wobbly: left endpoints of 6.07, 6.03, 6.03, 5.92, 6.03. Notice that not only do have these endpoints jumping around, there’s a mismatch between the one’s above and these.

bootstrap = do(100) * mean(resample(Data3350$Sleep))
qdata(~mean, p=c(0.025, 0.975), data=bootstrap)

    2.5%    97.5% 
6.030833 6.658182 

Starting to stabilize: 6.05, 6.06, 6.03, 6.04, 5.99. Except for the last one, we had a tight grouping.

bootstrap = do(500) * mean(resample(Data3350$Sleep))
qdata(~mean, p=c(0.025, 0.975), data=bootstrap)

    2.5%    97.5% 
5.984697 6.686742 

For the last one, our left endpoints were 6.04, 6.04, 6.03, 6.04, 6.05.

bootstrap = do(1000) * mean(resample(Data3350$Sleep))
qdata(~mean, p=c(0.025, 0.975), data=bootstrap)

    2.5%    97.5% 
6.048258 6.669697 

As we can see, a thousand trials is certainly enough, and the 50 and 100 were way too few. We did hit an outlier value with 500. Stats theory suggests about 500 is enough. However, there is no downside - as long as we have the computing power - to use a thousand or more trials. At a thousand plus, it appears we have stable endpoints for our confidence intervals.

Item 6.3

Use the VarsAth variable in Data3350 where Y / N responses indicate whether the participant’s is a varsity UNG athlete. Assuming the data frame is representative of the UNG Dahlonega campus, create a 95% confidence interval estimate for the percentage of students who are varsity athletes and interpret your findings. Hint: set success = “Y”.

We need to ensure the data pass verification: 10 or more successes and 10 or more failures in the sample.

tally(~VarsAth , data = Data3350)

VarsAth
  N   Y 
148  17 

Since we 17 successes and 148 failures, the data are appropriate for \(z\)-procedures. Even though RStudio is using the \(\chi^2\) statistic to evaluate it, we’re emulating a \(z\)-procedure, so I used the verification process for \(z\)-proportion procedures for my data checks.

Item 6.5

Use the CHS variable in Data3350 where numeric scores represent scores on the Coping Humor Scale. Assuming the data frame is representative of the UNG Dahlonega campus, create a 95% confidence interval estimate for the mean CHS score and interpret your findings.

We need to check the normality assumption.

favstats( ~CHS , data = Data3350)

With \(n = 163\) which is much greater than 40, the robustness of the \(t\)-statistic suggests good accuracy for the interval.

confint(t.test(~CHS, data = Data3350))

Item 6.6

Use the CHS variable in Data3350 to create a bootstrap confidence interval at the 95% level. Compare and contrast it with the results from the theoretical confidence interval. Use 500 resamplings.

Let’s create the bootstrapping resampling for the variable CHS and name the data frame boots.

boots = do(500) * mean(resample(Data3350$CHS))
qdata(~mean, p=c(0.025, 0.975), data=boots)

Notice there’s a good deal of agreement between the two different approaches to confidence intervals.

Item 7.2

Use the Data3350 data frame to build and evaluate a linear model for neuroticism (dependent variabele) vs. optimism (independent variable) using the Neuro and Opt variables. Be sure to check the linearity and normality assumptions and analyze all regression statistics. Construct a confidence interval for the slope of the regression line using an appropriate method.

xyplot(Neuro ~ Opt, data = Data3350)

We see and approximately downward trend, and certainly no non-linear pattern. Time to build the model and check the diagnostics.

modN = lm(Neuro ~ Opt, data = Data3350)
summary(modN)

Call:
lm(formula = Neuro ~ Opt, data = Data3350)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.153  -8.612   1.765   7.388  27.077 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  79.8572     4.2275  18.890  < 2e-16 ***
Opt          -1.8852     0.2096  -8.996 1.29e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.09 on 143 degrees of freedom
  (20 observations deleted due to missingness)
Multiple R-squared:  0.3614,    Adjusted R-squared:  0.357 
F-statistic: 80.94 on 1 and 143 DF,  p-value: 1.294e-15

With a multiple R-squared of 0.3614, we find that 36% of the variance in Neuro is accounted for by Opt. Since the absolute value of the correlation is \[r = \sqrt{R^2} = \sqrt{0.3614} = 0.601165534607565\] we see that \(r \approx -0.60\). Note that we have to manually add the negative sign for the correlation since the trend was downward sloping. Thus, we have moderate, negative correlation between Neuro and Opt.

Next, lets check the diagnostics.

qqmath(~resid(modN), type = c("p" , "r"))

While a few wiggles exist, the overall pattern lies reasonably close the 45 degree line and provides evidence for an approximately normal (bell-shaped) distribution of residuals. We can also check a fitted values plot, though they make more sense in the presence of multiple predictors.

plot(modN, which = 1 )

The red line is reasonably close to a straight line, and overall the model seems to work rather well and meet the assumptions for linear regression.

To analyze the slope coefficient of \(-1.89\) indicates that, for each 1 unit increase in Opt, we expect a \(1.89\) unit decrease in Neuro.

Item 7.4

Use the Perc data frame to create linear models for the 75th and 90th percentile wage earners from 2000 Q1 through 2019 Q4. Be sure to look carefully at all diagnostic plots. Create models for the Trump era and the last 3 years of the Obama era and conduct hypothesis tests that the Trump era growth was significantly greater than the historic trend as well as greater than the Obama era.

Full details available in my Wage Growth analysis. The summary chart below covers all 20 possible models.

US Weekly Wage Growth Per Quarter: All Wage Earners

\[\begin{array}{lccccc} &\textbf{10th Percentile}&\textbf{25th Percentile}&\textbf{Median}&\textbf{75th Percentile}&\textbf{90th Percentile}\\ \text{Historic} & \$2.03& \$2.77& \$4.05& \$7.49& \$12.02\\ \text{Trump Era} &\$4.95 & \$5.14 & \$6.90 & \$10.45 & \$20.72\\ \text{Obama Era} & \$1.49 & \$1.48 & \$3.20 & \$6.03 & \$8.34\\ \text{Obama Last 3 Yrs} & \$1.68 & \$2.06 & \$5.08 & \$6.40 & \$6.70 \end{array}\]

Item 7.5

The PercB data is similar to Perc in that it includes the 10th, 25th, 50th, 75th and 90th percentile wages, but PerB includes Black wage earners only. Use the PercB data frame to create a linear model for the median wage growth (50th percentile) from 2000 Q1 through 2019 Q4. Be sure to look carefully at all diagnostic plots. Create models for the Trump era and the last 3 years of the Obama era and conduct hypothesis tests that the Trump era growth was significantly greater than the historic trend as well as greater than the Obama era.

Full details available in my Wage Growth analysis. The summary chart below covers all 20 possible models.

US Weekly Wage Growth Per Quarter: African American Wage Earners Only

\[\begin{array}{lccccc} &\textbf{10th Percentile}&\textbf{25th Percentile}&\textbf{Median}&\textbf{75th Percentile}&\textbf{90th Percentile}\\ \text{Historic} & \$1.65 & \$2.25 & \$3.11 & \$5.26 & \$9.23\\ \text{Trump Era} & \$3.24 & \$4.66 & \$4.54 & \$9.82 & \$11.57\\ \text{Obama Era} & \$0.88 & \$1.32 & \$1.66 & \$3.79 & \$9.85\\ \text{Obama Last 3 Yrs} & \$1.68 & \$1.65 & \$0.92 & \$2.47 & \$4.60 \end{array}\]

Homework Set 4: Modules 8 and 9

Item 8.2

Using the HSAG variable from the Data3350 data frame, build a linear model for Aggressive Humor using the predictors Narcissism and Self-Defeating Humor Style: \[\text{HSAG} \sim \text{Narc} + \text{HSSD}\] Evaluate and analyze your model including all diagnostic plots.

The linear model:

modG = lm(HSAG ~ Narc + HSSD , data = Data3350)
summary(modG)

Call:
lm(formula = HSAG ~ Narc + HSSD, data = Data3350)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.5678  -3.3956  -0.5862   3.8085  15.6463 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 15.25834    2.00564   7.608 4.98e-12 ***
Narc         0.93980    0.18755   5.011 1.73e-06 ***
HSSD         0.31606    0.05922   5.337 4.07e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.926 on 130 degrees of freedom
  (32 observations deleted due to missingness)
Multiple R-squared:  0.296, Adjusted R-squared:  0.2851 
F-statistic: 27.32 on 2 and 130 DF,  p-value: 1.242e-10

Notice that R-squared = 0.296, so the two predictors together are accounting for 30% of the variance in Aggressive Humor.

qqmath(~resid(modG), type = c("p","r"))

plot(modG , which = 1)

The qq-plot is very much linear along the 45 degree line indicating the residuals follow an approximately normal distribution (normality assumption). The fitted values plot also shows the red fit-line lies very close to the dotted horizontal line indicating perfect linearity (linearity assumption). Thus, these data are well-suited for linear regression modeling.

Consider the interpretation of the model: Aggressive humor depends on another humor style – not surprising two types of humor are correlated – and also on Narcissism. That makes sense, too. Aggressive humor is comprised of sarcasm and put-downs, and narcissists would find it less irksome to belittle others.

Item 8.3

Add the Eating Attitudes variable Eat as the third predictor in your model for HSAG:\[\text{HSAG} \sim \text{Narc} + \text{HSSD} + \text{Eat}\] Evaluate and analyze your model including all diagnostic plots, and compare your three-predictor model with your results from the two-predictor model. Is the correlation between Eat and HSAG positive or negative? How can you tell from the model summary statistics output? If higher scores on the Eat variable indicate higher levels of being calorie conscious, knowing the fat content of food items and thinking about burning calories when working out, does this relationship between Eat and HSAG make sense?

The model.

modG2 = lm(HSAG ~ Narc + HSSD + Eat, data = Data3350)
summary(modG2)

Call:
lm(formula = HSAG ~ Narc + HSSD + Eat, data = Data3350)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.3174  -3.0466  -0.3175   3.6231  12.8209 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 19.63898    2.54922   7.704 3.28e-12 ***
Narc         1.00601    0.18748   5.366 3.69e-07 ***
HSSD         0.30799    0.05921   5.201 7.70e-07 ***
Eat         -0.15245    0.05825  -2.617  0.00994 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.828 on 127 degrees of freedom
  (34 observations deleted due to missingness)
Multiple R-squared:  0.3196,    Adjusted R-squared:  0.3035 
F-statistic: 19.89 on 3 and 127 DF,  p-value: 1.248e-10

We have increased R-squared to 0.3196, so we are now accounting for an additional 2% of the variance in the Aggressive Humor variable. Note the \(p\)-values for the significance of the correlation. While the Eating Attitudes variable is much less significant than the other two, it’s still significant at the 0.01 level which is why it’s able to contribute additional R-squared to the overall model.

qqmath(~resid(modG2), type = c("p", "r"))

plot(modG2 , which = 1)

The linearity assumption checks out – the fitted values plot shows nearly perfect linearity. The qq-plot wobbles a bit, but the residuals do seem to be approximately normally distributed. Our three-predictor model is excellent, and one would very much trust research conclusions based on a model with summary statistics and diagnostics like this one.

The addition of the Eating Attitudes variable is interesting. The questions are about knowing the number of calories in foods we eat, thinking about how many calories are burned when exercising and other health conscious behaviors. Notice this coefficient is negatively correlated with Aggressive Humor. Those who are more health conscious may – I’m totally guessing here – have an aversion to humor based on insults and put-downs because they are more conscious about weight and body image. Perhaps their empathy makes them resistant to aggressive humor. Whatever the reasons, it is an interesting relationship among personality variables and bears consideration.

Item 9.2

A supposedly fair die is rolled 50 times with the following results: 9 ones, 15 twos, 9 threes, 8 fours, 6 fives and 13 sixes. Test at the 0.05 level whether the die is actually fair using \(\chi^2\) GOF.

I find it hilarious that this fails:

observed = c(9 , 15 , 9 , 8 , 6 , 13)
expected = c(0.1667 , 0.1667 , 0.1667 , 0.1667 , 0.1667 , 0.1667 )
xchisq.test(x = observed,
           p = expected)

Error in chisq.test(x = x, y = y, correct = correct, p = p, rescale.p = rescale.p,  : 
  probabilities must sum to 1.

What RStudio is expecting is that we will leave the expected parameter to its default setting when each probability is equally likely. So this works:

xchisq.test(x = observed)

    Chi-squared test for given probabilities

data:  x
X-squared = 5.6, df = 5, p-value = 0.3471

  9.00    15.00     9.00     8.00     6.00    13.00  
(10.00)  (10.00)  (10.00)  (10.00)  (10.00)  (10.00) 
 [0.10]   [2.50]   [0.10]   [0.40]   [1.60]   [0.90] 
<-0.32>  < 1.58>  <-0.32>  <-0.63>  <-1.26>  < 0.95> 
           
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

Also, if you let RStudio create the rounding error, the values will (in R’s complicated mind) sum to 1. Let’s create a variable for the probability of a dice roll:

pDice = 1/6
pDice

[1] 0.1666667

observed = c(9 , 15 , 9 , 8 , 6 , 13)
expected = c(pDice , pDice , pDice , pDice , pDice , pDice )
xchisq.test(x = observed,
           p = expected)

    Chi-squared test for given probabilities

data:  x
X-squared = 5.6, df = 5, p-value = 0.3471

  9.00    15.00     9.00     8.00     6.00    13.00  
(10.00)  (10.00)  (10.00)  (10.00)  (10.00)  (10.00) 
 [0.10]   [2.50]   [0.10]   [0.40]   [1.60]   [0.90] 
<-0.32>  < 1.58>  <-0.32>  <-0.63>  <-1.26>  < 0.95> 
           
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

So we’ve tested the probability model (null hypothesis) using \(\chi^2\) GOF and found a \(p\)-value of 0.3471. The expected cell counts are all equal to 10 which is greater than 5, so the data pass the verification check. The statistical conclusion? We fail to reject the null. We have no evidence the probability model is incorrect.