<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/"></div></div><div><a data-url="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/" data-counturl="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/" data-count="horizontal" data-text="Raccoon | Ch 2.2 – 2 sample t-test and paired t" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.quantide.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/"> R blog | Quantide - R training &amp; consulting</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><em>This article is part of Quantide’s web book “<a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-statistical-models-with-r-web-book/">Raccoon – Statistical Models with R</a>“. Raccoon is Quantide’s third web book after “<a target="_blank" rel="nofollow" href="http://www.quantide.com/rabbit-introduction-to-r/">Rabbit – Introduction to R</a>” and “<a target="_blank" rel="nofollow" href="http://www.quantide.com/ramarro-r-for-developers/">Ramarro – R for Developers</a>“. See the <a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-statistical-models-r-free-web-book/">full project here</a>.</em></p>
<p><em>The second chapter of Racooon is focused on T-test and Anova. Through example it shows theory and R code of:</em></p>
<ul>
<li><a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-1-1-sample-t-test-with-r/">1-sample t-test</a></li>
<li><a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/">2-sample t-test and Paired t</a></li>
<li>1-way Anova</li>
<li>2-way Anova</li>
<li>Unbalanced Anova with one obs dropped</li>
<li>Fixed effects Nested Anova</li>
</ul>
<p><em>This post is the second section of the chapter, about 2-sample t-test and paired t. </em></p>
<p><em>Throughout the web-book we will widely use the <strong>package qdata</strong>, containing about 80 datasets. You may find it here: <a target="_blank" rel="nofollow" href="https://github.com/quantide/qdata">https://github.com/quantide/qdata</a>.</em></p>
<pre>require(nortest)
require(car)
require(ggplot2)
require(dplyr)
require(gridExtra)
require(tidyr)
require(qdata)</pre>

<h2>Example: Chicken hormones (2 -sample \(t\)-test)</h2>
<h3>Data description</h3>
<p>To compare two growth hormones, 18 chicken were tested: 10 being assigned to hormone A and 8 to hormone B. The gains in weights (grams) over the period of experiment were measured.</p>
<h3>Data loading</h3>
<p>Let us load and have a look at the data</p>
<pre>data(hormones)
head(hormones)</pre>

<pre>## # A tibble: 6 × 2
##   hormone  gain
##     
## 1       A   615
## 2       A   645
## 3       A   840
## 4       A   645
## 5       A   365
## 6       A   450</pre>

<pre>str(hormones)</pre>

<pre>## Classes 'tbl_df', 'tbl' and 'data.frame':    18 obs. of  2 variables:
##  $ hormone: Factor w/ 2 levels &quot;A&quot;,&quot;B&quot;: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gain   : int  615 645 840 645 365 450 530 756 785 560 ...</pre>

<h2>Descriptives</h2>
<p>We begin by calculating some descriptive statistics:</p>
<pre>hormones %>% group_by(hormone) %>%
  summarise(min = min(gain), 
    first_qu = quantile(gain, 0.25),
    median = median(gain),
    mean = mean(gain),
    third_qu = quantile(gain, 0.75),
    max = max(gain),
    sd = sd(gain))</pre>

<pre>## # A tibble: 2 × 8
##   hormone   min first_qu median    mean third_qu   max       sd
##                       
## 1       A   365    537.5  630.0 619.100   728.25   840 149.4482
## 2       B   560    742.0  827.5 822.625   916.25  1035 154.5574</pre>
<p>And now we draw a boxplot and a stripchart (in this case in vertical) of the data:</p>
<pre>ggp &lt;- ggplot(data = hormones, mapping = aes(x=hormone, y=gain, fill=hormone)) +
  geom_boxplot()

print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova5" class="wp-image-2891" src="https://i1.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova5-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Box-and-whiskers plot of gain by hormone</p>

<pre>hormones_mean &lt;- hormones %>% group_by(hormone) %>% summarise(gain=mean(gain))

ggp &lt;- ggplot(data = hormones, mapping = aes(x=hormone, y=gain)) +
  geom_jitter(position = &quot;dodge&quot;, color=&quot;darkblue&quot;) +
  geom_point(data=hormones_mean, mapping = aes(x=hormone,y=gain), colour=&quot;red&quot;, group=1) +
  geom_line(data=hormones_mean, mapping = aes(x=hormone,y=gain), colour=&quot;red&quot;, group=1) 
  
print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova6" class="wp-image-2892" src="https://i2.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova6-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Stripchart with connect line of gain by hormone</p>
<p>The <code>geom_point()</code> function and <code>geom_line()</code> functions add the red points of two group averages to second graph and then link them with a segment.<br>
<code>group_by()</code> and <code>summarise()</code> functions applies separately the <code>mean()</code> function to values of <code>gain</code> corresponding to <code>hormone</code> variable levels. In this case, if the segment is near to parallel to the x-axis, then the two means are probably equal. Anyway, both graphs and statistics show that hormone B tends to generate heavier chickens.</p>
<h3>Inference and models</h3>
<p>Let us check the normality of <code>gain</code> in the two groups through q-q plots</p>
<pre>hormones_a &lt;- hormones %>% filter(hormone ==&quot;A&quot;) 
hormones_b &lt;- hormones %>% filter(hormone ==&quot;B&quot;)

ggp1 &lt;- ggplot(data = hormones_a, mapping = aes(sample = gain)) + 
  stat_qq(color=&quot;#F8766D&quot;, size=2) +
  geom_abline(intercept=mean(hormones_a$gain),slope=sd(hormones_a$gain), color=&quot;red&quot;, linetype=2) +
  ylab(&quot;hormone A&quot;) + ggtitle(&quot;q-q plot of hormone A&quot;)

ggp2 &lt;- ggplot(data = hormones_b, mapping = aes(sample = gain)) + 
  stat_qq(color=&quot;#00C094&quot;, size=2) +
  geom_abline(mapping = aes(intercept=mean(hormones_b$gain),slope=sd(hormones_b$gain)), color=&quot;red&quot;, linetype=2) +
  ylab(&quot;hormone B&quot;) + ggtitle(&quot;q-q plot of hormone B&quot;)

grid.arrange(ggp1, ggp2, nrow = 2)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova7" class="wp-image-2893" src="https://i2.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova7-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Normal probability plot of gain by hormone</p>
<p>From graph, no evidence of non-normality is evident. Now let us perform Anderson-Darling’s test for <code>gain</code> in the two groups, exploiting the function <code>tapply</code>:</p>
<pre>tapply(hormones$gain, hormones$hormone, ad.test)</pre>

<pre>## $A
## 
##  Anderson-Darling normality test
## 
## data:  X[[i]]
## A = 0.15886, p-value = 0.926
## 
## 
## $B
## 
##  Anderson-Darling normality test
## 
## data:  X[[i]]
## A = 0.13781, p-value = 0.9557</pre>
<p>The last line of above code, returns a 2-element list, with AD test for each sub-group. There is no evidence to reject the normality hypothesis within the two groups.<br>
The next line of code tests the homoscedasticity assumption, that is if the variances within the two groups are equal:</p>
<pre>var.test(gain ~ hormone, data = hormones)</pre>

<pre>## 
##  F test to compare two variances
## 
## data:  gain by hormone
## F = 0.93498, num df = 9, denom df = 7, p-value = 0.9031
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1938497 3.9241513
## sample estimates:
## ratio of variances 
##          0.9349792</pre>

<pre>bartlett.test(gain ~ hormone,data=hormones)</pre>

<pre>## 
##  Bartlett test of homogeneity of variances
## 
## data:  gain by hormone
## Bartlett's K-squared = 0.0083868, df = 1, p-value = 0.927</pre>

<pre>leveneTest(gain ~ hormone,data=hormones)</pre>

<pre>## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0321 0.8601
##       16</pre>
<p><span>All three previous tests accept \(H_0\), i.e., there is no evidence to say that the variability changes across levels.</span></p>
<p><span>Now the \(t\)-test compares the two group means.  </span></p>
<p><span>The null hypothesis is \(H_0: beta_1 = beta_2\), where \(beta_1\) and \(beta_2\) are the expectations of dependent variable for hormones A and B, respectively.</span></p>
<pre>t.test(gain ~ hormone, data = hormones, var.equal = TRUE)</pre>

<pre>## 
##  Two Sample t-test
## 
## data:  gain by hormone
## t = -2.8283, df = 16, p-value = 0.01211
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -356.07303  -50.97697
## sample estimates:
## mean in group A mean in group B 
##         619.100         822.625</pre>
<p>The \(t\)-test rejects the null hypothesis.</p>
<p>Now the same test may be conducted using a model approach.<br>
The following lines of code perform the model fitting and print the main results.</p>
<pre>mod &lt;- lm(gain ~ hormone, data = hormones) 
mod</pre>

<pre>##</pre>

<pre>## Call:
## lm(formula = gain ~ hormone, data = hormones)
## 
## Coefficients:
## (Intercept)     hormoneB  
##       619.1        203.5</pre>
<p>The parameters shown are the model intercept and the “slope” coefficient for hormone <code>B</code>. The intercept in this case is the estimate of mean growth for <code>A</code> hormone, whereas the slope in this case is the estimate of difference between the mean growth obtained with hormon <code>B</code> and the mean growth obtained with hormone <code>A</code>.<br>
In other words, the mean growth estimate for <code>hormoneA</code> is <code>619.1</code>, while the mean growth estimate for <code>hormoneB</code> is <code>619.1 + 203.5 = 822.6</code>. We will return to parameterization issue later.</p>
<p>Using <code>summary</code> the significancy of parameters is also shown:</p>
<pre>summary(mod)</pre>

<pre>## 
## Call:
## lm(formula = gain ~ hormone, data = hormones)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -262.62  -83.48   10.90  120.77  220.90 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   619.10      47.97  12.905 7.12e-10 ***
## hormoneB      203.53      71.96   2.828   0.0121 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 151.7 on 16 degrees of freedom
## Multiple R-squared:  0.3333, Adjusted R-squared:  0.2916 
## F-statistic: 7.999 on 1 and 16 DF,  p-value: 0.01211</pre>
<p>The test results are equal to the ones of \(t\)-test, and the \(t\)-test value is equal to minus \(t\)-test value of <code>hormoneB</code> coefficient.</p>
<p>Why these results? “Behind the scene”, the model has been built using the following **model matrix** (\(underline{X}\))</p>
<pre>model.matrix(mod)</pre>

<pre>##    (Intercept) hormoneB
## 1            1        0
## 2            1        0
## 3            1        0
## 4            1        0
## 5            1        0
## 6            1        0
## 7            1        0
## 8            1        0
## 9            1        0
## 10           1        0
## 11           1        1
## 12           1        1
## 13           1        1
## 14           1        1
## 15           1        1
## 16           1        1
## 17           1        1
## 18           1        1
## attr(,&quot;assign&quot;)
## [1] 0 1
## attr(,&quot;contrasts&quot;)
## attr(,&quot;contrasts&quot;)$hormone
## [1] &quot;contr.treatment&quot;</pre>
<p>What this matrix means?</p>
<p>For modeling purposes, the independent variable <code>hormone</code> has been “split” in two dummy variables: <code>hormoneA</code> and <code>hormoneB</code>:</p>
<ul>
<li><code>hormoneA</code> is <code>1</code> when data contain measures relative to hormone <code>A</code>, and is <code>0</code> when data contain measures relative to hormone <code>B</code></li>
<li><code>hormoneB</code> is <code>1</code> when data contain measures relative to hormone <code>B</code>, and is <code>0</code> when data contain measures relative to hormone <code>A</code></li>
</ul>
<p>Obviously, when <code>hormoneA</code> is <code>0</code>, <code>hormoneB</code> is <code>1</code> and vice versa. Then, one dummy variable (in this case, <code>hormoneA</code>), can be removed from the model, because it does not add any useful information, and makes the model mathematically non-manageable. Consequently, the model matrix will contains only the “all 1s” intercept column (that refers to the case when hormone <code>A</code> is used) and the “displacement” <code>hormoneB</code>column which represents the constant added to mean growth when using hormone <code>B</code> instead of hormone <code>A</code>.</p>
<p>In other words, the model actually estimated is the following:</p>
<p><span>\( growth_i = beta_0 + beta_1 cdot delta_B(hormone_i) + varepsilon_i \)</span></p>
<p><span>Where \(delta_B(hormone_i)\) represents the function that returns 1 only when the hormone of \(i\)-th unit is <code>B</code>, otherwise, returns 0.  </span></p>
<p><span>\(beta_0\) in this manner is the value of mean growth when using hormone <code>A</code> (i.e., when \(delta_B(hormone_i)==0\)), while \(beta_1\) is the “jump” in mean of growth when using hormone <code>B</code> instead of hormone <code>A</code>.</span></p>
<p>Another approach available to manage such type of test is through the analysis of variance model, i.e. by using the <code>aov()</code> function, that directly performs one analysis of variance:</p>
<pre>fm &lt;- aov(mod)</pre>
<p>Printing a summary method for analysis of variance objects reflects the different approach to what is essentially the same information as returned by <code>lm()</code>. While a printed linear-model object shows individual coefficients, printed <code>aov</code> objects show terms, which may correspond to several coefficients. However, methods for linear models can always be applied explicitly on <code>aov</code> objects, such as those for the coefficients or residuals.</p>
<pre>summary(fm)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## hormone      1 184100  184100   7.999 0.0121 *
## Residuals   16 368229   23014                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>
<p>Above output shows that <code>hormone</code>, with all its levels, is globally significant (p=0.012) in explaining the dependent variable variability.</p>
<p>In this case, since the independent variable (factor) contains two levels only, the resulting p-value for <code>hormone</code> obtained from <code>aov</code> is identical to the p-value of hormone <code>B</code> obtained using <code>lm()</code> (actually, in this example they are the same test).</p>
<p>If the independent variable (factor) contains more than two levels, then <code>aov()</code> returns the global significancy (not the “by level significancy”) of independent variable effect on dependent variable.</p>
<p>A few more concise code to obtain the same results as above is</p>
<pre>fm &lt;- aov(gain ~ hormone, data = hormones) 
summary(fm)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## hormone      1 184100  184100   7.999 0.0121 *
## Residuals   16 368229   23014                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>
<p>Next lines of code show that actually <code>aov</code> objects are <code>lm</code> objects with also analysis of variance results added.</p>
<pre>class(mod)</pre>

<pre>## [1] &quot;lm&quot;</pre>

<pre>class(fm)</pre>

<pre>## [1] &quot;aov&quot; &quot;lm&quot;</pre>
<p><span>The following commands are useful to read/define the type of contrasts (i.e., the type of coding of regressors matrix \(underline{X}\)) used when unordered or ordered factors are used.</span></p>
<pre>options(&quot;contrasts&quot;)</pre>

<pre>## $contrasts
##         unordered           ordered 
## &quot;contr.treatment&quot;      &quot;contr.poly&quot;</pre>

<pre>options(contrasts = c(&quot;contr.treatment&quot;, &quot;contr.poly&quot;))</pre>
<p>The first line of code simply shows the current predefined type of contrasts for unordered and ordered factors, while the second line of code sets the predefined contrasts.</p>
<p>To extract coefficients from an already calculated (<code>lm</code> or <code>aov</code>) model, one must use the <code>coefficients</code> method:</p>
<pre>coefficients(fm)</pre>

<pre>## (Intercept)    hormoneB 
##     619.100     203.525</pre>
<p>Next line of code shows that the <code>hormoneB</code> coefficient is simply the difference between average growths for <code>hormoneA</code> and <code>hormoneB</code> chickens.</p>
<pre>diff(tapply(hormones$gain, hormones$hormone, mean))</pre>

<pre>##       B 
## 203.525</pre>
<p>R can use several types of contrasts: <code>contr.treatment</code> sets to zero the coefficient of first factor level, that becomes the intercept, and consequently all the means of other levels are relative to the first level. Another option is <code>contr.sum</code>, that sets to zero the sum of coefficients, and then the reference value of coefficients becomes the gran mean, calculated on all factor levels. R uses <code>contr.treatment</code> as default.</p>
<p>Let’s try to change type of contrasts:</p>
<pre>options(contrasts = c(&quot;contr.sum&quot;, &quot;contr.poly&quot;))</pre>

<pre>fm &lt;- aov(gain~hormone, data = hormones)</pre>
<p>The <code>aov()</code> results are the same:</p>
<pre>summary(fm)</pre>

<pre>##             Df Sum Sq Mean Sq F value Pr(>F)  
## hormone      1 184100  184100   7.999 0.0121 *
## Residuals   16 368229   23014                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>
<p>but, on the whole, the results are changed:</p>
<pre>coefficients(fm)</pre>

<pre>## (Intercept)    hormone1 
##    720.8625   -101.7625</pre>
<p><span>The values of the coefficients are not the same because with <code>contr.treatment</code> the intercept represents the mean in the first group (hormone A) and the “slope” represents the difference between the mean in the second group (hormone B) and the mean in the first group; conversely, with <code>contr.sum</code> the intercept represents the grand mean, and the slope represents the difference between the mean in the first group and the grand mean.<br>
The differences in coefficients are perhaps more perceptible looking at the model matrix:</span></p>
<pre>model.matrix(fm)</pre>

<pre>##    (Intercept) hormone1
## 1            1        1
## 2            1        1
## 3            1        1
## 4            1        1
## 5            1        1
## 6            1        1
## 7            1        1
## 8            1        1
## 9            1        1
## 10           1        1
## 11           1       -1
## 12           1       -1
## 13           1       -1
## 14           1       -1
## 15           1       -1
## 16           1       -1
## 17           1       -1
## 18           1       -1
## attr(,&quot;assign&quot;)
## [1] 0 1
## attr(,&quot;contrasts&quot;)
## attr(,&quot;contrasts&quot;)$hormone
## [1] &quot;contr.sum&quot;</pre>
<p>Now, firstly we will restore predefined contrasts</p>
<pre>options(contrasts = c(&quot;contr.treatment&quot;, &quot;contr.poly&quot;))</pre>
<p>and then we will fit a model without intercept to data. In general, removing the intercept is a discouraged practice, but sometimes no intercept models can help to enhance coefficients meaning:</p>
<pre>mod &lt;- lm(gain ~ hormone-1, data = hormones)
coefficients(mod)</pre>

<pre>## hormoneA hormoneB 
##  619.100  822.625</pre>
<p>The above coefficients represent, respectively, the mean in the group <code>hormoneA</code> and in the group <code>hormoneB</code>.</p>
<p>The model matrix can be useful to understand the meaning of the coefficients:</p>
<pre>model.matrix(mod)</pre>

<pre>##    hormoneA hormoneB
## 1         1        0
## 2         1        0
## 3         1        0
## 4         1        0
## 5         1        0
## 6         1        0
## 7         1        0
## 8         1        0
## 9         1        0
## 10        1        0
## 11        0        1
## 12        0        1
## 13        0        1
## 14        0        1
## 15        0        1
## 16        0        1
## 17        0        1
## 18        0        1
## attr(,&quot;assign&quot;)
## [1] 1 1
## attr(,&quot;contrasts&quot;)
## attr(,&quot;contrasts&quot;)$hormone
## [1] &quot;contr.treatment&quot;</pre>
<p><span>Next code shows what can happen when one uses a no-intercept model and does not address for that in analysis of variance:</span></p>
<pre>fm &lt;- aov(mod) 
summary(fm)</pre>

<pre>##           Df  Sum Sq Mean Sq F value   Pr(>F)    
## hormone    2 9246543 4623272   200.9 4.63e-12 ***
## Residuals 16  368229   23014                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<pre>coefficients(fm)</pre>

<pre>## hormoneA hormoneB 
##  619.100  822.625</pre>
<p>The significancy of factor is much more strong than in previous calculations, but this is wrong, because this result actually compares the model with <code>0</code> intercept with the model that uses one mean for <code>hormoneA</code>and one mean for <code>hormoneB</code>.</p>
<p>Using the following code gives results equivalent to t-test:</p>
<pre>anova(mod,lm(gain ~ 1, data = hormones))</pre>

<pre>## Analysis of Variance Table
## 
## Model 1: gain ~ hormone - 1
## Model 2: gain ~ 1
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     16 368229                              
## 2     17 552328 -1   -184100 7.9994 0.01211 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>
<p><code>anova()</code> performs an analysis of variance to compare two nested models (see the paragraph on <em>Comparing two nested models</em> in <em>Some Theory on Linear Models</em> chapter).<br>
In this case we are comparing a model with the group variable <code>hormone</code> against the model with grand mean only.<br>
The result is the same as seen individually in <code>lm()</code> and <code>aov()</code> outputs.</p>
<h3>Model diagnostics</h3>
<p>Since the assumpions on linear models, the residuals should arise from a gaussian distribution; then, a Normal Probability Plot of residuals will be plotted to check this assumption.</p>
<pre>res &lt;- cbind(hormones,res=fm$residuals)

ad.test(res$res)</pre>

<pre>## 
##  Anderson-Darling normality test
## 
## data:  res$res
## A = 0.18963, p-value = 0.8859</pre>

<pre>ggp &lt;- ggplot(data = res, mapping = aes(sample = res, color=hormone)) + 
  stat_qq(size=2) +
  geom_abline(mapping = aes(intercept=mean(res), slope=sd(res)), colour=&quot;red&quot;, linetype=2)

print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova8" class="wp-image-2894" src="https://i0.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova8-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Overlaid Normal probability plots of residuals</p>
<p>The two colors correspond to two groups.<br>
Residuals seem gaussian, irrespective of group.</p>
<p>The Anderson-Darling (AD) test confirms that:</p>
<pre>ad.test(fm$residuals)</pre>

<pre>## 
##  Anderson-Darling normality test
## 
## data:  fm$residuals
## A = 0.18963, p-value = 0.8859</pre>
<p>This confirms that the analysis and its results may be considered correct.</p>
<h2>Example: Park assistant <span>(paired \(t\))</span></h2>
<h3>Data description</h3>
<p>In an experiment to investigate if an optional equipment of cars may help to park the cars themselves, 20 drivers were asked to park the car with and without the equipment. The parking time (secs) for each driver and equipment is recorded.</p>
<h3>Data loading</h3>

<pre>data(carctl)
str(carctl)</pre>

<pre>## Classes 'tbl_df', 'tbl' and 'data.frame':    20 obs. of  2 variables:
##  $ Car_A: num  35.4 28.9 32.1 40.6 18.3 ...
##  $ Car_B: num  28.6 24.4 37.3 39.9 21.7 ...</pre>

<h3>Descriptives</h3>
<p>A simple boxplot:</p>
<pre>carctl2 % gather(car, value, Car_A, Car_B)</pre>

<pre>ggp &lt;- ggplot(data = carctl2, mapping = aes(x=car, y=value, fill=car)) +
  geom_boxplot()

print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova9" class="wp-image-2895" src="https://i0.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova9-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Box-and-whiskers plot of parking times for the two cars</p>
<p><span>A stripchart with mean points and mean line:</span></p>
<pre>carctl_mean % group_by(car) %>% summarise(value=mean(value))</pre>

<pre>ggp &lt;- ggplot(data = carctl2, mapping = aes(x=car, y=value)) +
  geom_jitter(position = &quot;dodge&quot;, color=&quot;darkblue&quot;) +
  geom_point(data=carctl_mean, mapping = aes(x=car,y=value), colour=&quot;red&quot;, group=1) +
  geom_line(data=carctl_mean, mapping = aes(x=car,y=value), colour=&quot;red&quot;, group=1) 
  
print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova10" class="wp-image-2896" src="https://i1.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova10-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Stripchart with connect line of parking times for the two cars</p>
<p>And now the \(t\)-test.</p>
<p>Before testing the means, a check for normality is made by using qq-plot and Anderson-Darling test</p>
<pre>ggp &lt;- ggplot(data = carctl2, mapping = aes(sample = value, color=car)) + 
  stat_qq() +
  geom_abline(data= carctl2 %>% filter(car==&quot;Car_A&quot;), mapping=aes(intercept=mean(value), 
    slope=sd(value), color=car), size=1) +
  geom_abline(data= carctl2 %>% filter(car==&quot;Car_B&quot;), mapping=aes(intercept=mean(value),
    slope=sd(value), color=car), size=1)</pre>

<pre>print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova11" class="wp-image-2897" src="https://i1.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova11-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Normal probability plot of parking times of two cars</p>
<p>There is not evidence of non-normality within groups, and the AD tests confirm above results:</p>
<pre>lapply(X=carctl, FUN=ad.test)</pre>

<pre>## $Car_A
## 
##  Anderson-Darling normality test
## 
## data:  X[[i]]
## A = 0.19988, p-value = 0.8647
## 
## 
## $Car_B
## 
##  Anderson-Darling normality test
## 
## data:  X[[i]]
## A = 0.36913, p-value = 0.3928</pre>
<p>and finally the two independent samples \(t\)-test:</p>
<pre>t.test(x=carctl$Car_A, y=carctl$Car_B)</pre>

<pre>## 
##  Welch Two Sample t-test
## 
## data:  carctl$Car_A and carctl$Car_B
## t = 1.2205, df = 36.778, p-value = 0.2301
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.056978  8.285186
## sample estimates:
## mean of x mean of y 
##  34.56418  31.45008</pre>
<p>Following this result, the null hypothesis is not rejected (at a 5% level).</p>
<p>But if we produce a scatterplot of <code>Car_A</code> Vs <code>Car_B</code>,</p>
<pre>ggp &lt;- ggplot(data = carctl, mapping = aes(x=Car_A, y=Car_B)) + 
  geom_point(col=&quot;darkblue&quot;)

print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova12" class="wp-image-2898" src="https://i0.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova12-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Scattrplot of parking times of Car_A Vs. Car_B</p>
<p>we can see that data of two cars are clearly NOT independent!!<br>
Measures on <code>Car_A</code> are strongly related to measures on <code>Car_B</code>, because each driver drove either cars.</p>
<p>Then the correct test for such as data is:</p>
<pre>carctl % mutate(car_diff = Car_A - Car_B)
t.test(x=carctl$car_diff, mu=0)</pre>

<pre>## 
##  One Sample t-test
## 
## data:  carctl$car_diff
## t = 3.8735, df = 19, p-value = 0.001023
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.431407 4.796801
## sample estimates:
## mean of x 
##  3.114104</pre>
<p>or, equivalently, a paired \(t\)-test(x=carctl\(Car_A, y=carctl\)Car_B, paired=TRUE)</p>
<pre>t.test(x=carctl$Car_A, y=carctl$Car_B, paired=TRUE)</pre>

<pre>## 
##  Paired t-test
## 
## data:  carctl$Car_A and carctl$Car_B
## t = 3.8735, df = 19, p-value = 0.001023
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.431407 4.796801
## sample estimates:
## mean of the differences 
##                3.114104</pre>
<p>In either cases the results are the same: a difference statistically significant exists between mean parking times of two cars.</p>
<p>In this case the Normality test must be performed on differences between measures within same driver.</p>
<p>Check of normality by Normal probability plot:</p>
<pre>ggp &lt;- ggplot(data = carctl, mapping = aes(sample = car_diff)) + 
  stat_qq(color=&quot;darkblue&quot;, size=2) +
  geom_abline(mapping = aes(intercept=mean(car_diff), slope=sd(car_diff)), colour=&quot;red&quot;, linetype=2)

print(ggp)</pre>
<p><img sizes="(max-width: 450px) 100vw, 450px" alt="anova-13" class="wp-image-2899" src="https://i2.wp.com/www.quantide.com/wp-content/uploads/2016/09/Anova-13-300x300.png?resize=450%2C450" data-recalc-dims="1"></p>
<p>Normal probability plot of difference of parking times of two cars</p>
<p>And the Anderson-Darling test:</p>
<pre>ad.test(x=carctl$car_diff)</pre>

<pre>## 
##  Anderson-Darling normality test
## 
## data:  carctl$car_diff
## A = 0.29177, p-value = 0.5704</pre>
<p>For either tests, the Normality assumption is not rejected.</p>
<p>Now we try the model approach to the same test. The model approach in this case follows the same procedure seen for the one-sample t-test, but using the differences between times as dependent variable.</p>
<p>Fitting of model:</p>
<pre>mod &lt;- lm(car_diff ~ 1, data = carctl) 
mod</pre>

<pre>## 
## Call:
## lm(formula = car_diff ~ 1, data = carctl)
## 
## Coefficients:
## (Intercept)  
##       3.114</pre>

<pre>summary(mod)</pre>

<pre>## 
## Call:
## lm(formula = car_diff ~ 1, data = carctl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2986 -2.1332  0.4141  2.3877  5.6064 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.114      0.804   3.873  0.00102 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.595 on 19 degrees of freedom</pre>
<p>or, equivalently one may use also:</p>
<pre>mod &lt;- lm(I(Car_A-Car_B)~1, data=carctl)
summary(mod)</pre>

<pre>## 
## Call:
## lm(formula = I(Car_A - Car_B) ~ 1, data = carctl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2986 -2.1332  0.4141  2.3877  5.6064 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.114      0.804   3.873  0.00102 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.595 on 19 degrees of freedom</pre>
<p>q-q plot (or Normal probability plot) of residuals shall be identical to previous one.</p>
<p>The Anderson-Darling test too is identical to previuosly made test:</p>
<pre>ad.test(mod$residuals)</pre>

<pre>## 
##  Anderson-Darling normality test
## 
## data:  mod$residuals
## A = 0.29177, p-value = 0.5704</pre>
<p>equivalently:</p>
<pre>ad.test(residuals(mod))</pre>

<pre>## 
##  Anderson-Darling normality test
## 
## data:  residuals(mod)
## A = 0.29177, p-value = 0.5704</pre>
<p>Note that for repeated measures experiments, several analytical options are also available: for example through the <code>Error()</code> function in formula specification, or using mixed-model methods (<code>lme4</code> or <code>nlme</code>packages are the most used and known).<br>
In this manual we won’t deep these arguments, because they would require a distinct section for each of them.</p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.quantide.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/">Raccoon | Ch 2.2 – 2 sample t-test and paired t</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.quantide.com/">Quantide – R training &amp; consulting</a>.</p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/"></div></div><div><a data-url="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/" data-counturl="https://www.r-bloggers.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/" data-count="vertical" data-text="Raccoon | Ch 2.2 – 2 sample t-test and paired t" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.quantide.com/raccoon-ch-2-2-2-sample-t-test-and-paired-t/"> R blog | Quantide - R training &amp; consulting</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/zVGP3a4DdW4"></a></div>