In ehassler/MontgomeryDAE: MontgomeryDAE - Data and examples from Montgomery's "Design and Analysis of Experiments"
library(knitr)
library(kableExtra)
library(MontgomeryDAE)
source('_format.R')
Chapter 2 begins with several methods of visualization that are easily reproducible within R. These visualizations are based on the data in Table 2.1 which records the results of an experiment where a researcher was comparing the strength of a new mortar formulation to that of the original unmodified formulation.
kable(Table2.1) %>% kable_styling()
In R, the object Table2.1 is a data.frame that has a column ModifiedMortar and a column UnmodifiedMortar. The structure of this object is given below.
str(Table2.1)
\noindent We can use the R command attach to take the columns of the data.frame and make them accessible by name (e.g. we can use ModifiedMortar instead of Table2.1$ModifiedMortar):
attach(Table2.1)
all(ModifiedMortar == Table2.1$ModifiedMortar)
\noindent When we're all done with the data we should call detach(Table2.1).
With the data loaded let's move on to some of the visualizations in Chapter 2.
Visualizations
We begin with the dot plot of Figure 2.1:
stripchart(
list(
'Modified'=ModifiedMortar,
'Unmodified'=UnmodifiedMortar
),
main="Dot Diagram for Table 2.1",
xlab='Strength'
)
\noindent The graph seems to suggest that the modified formulation has a lower average strength than that of the unmodified formulation.
Figure 2.2 in the text is a histogram of 200 observations of metal recovery for which no data is given. Though we lack the full data set, let's construct a histogram of the 20 data points from combining the modified and unmodified mortar formulation data:
data <- c(ModifiedMortar, UnmodifiedMortar)
hist(data)
\noindent Though it was easy to construct, the histogram is not giving us a lot of information in this case due to the fact that we have such a small amount of data.
The dot plot of Figure 2.1 can become very difficult to interpret when the number of observations becomes large. In these cases we can use a box plot which is constructed based on summary statistics of the data. The box plots of Figure 2.3 are presented next.
bp.returned <- boxplot(
Table2.1,
main="Box plots for Table 2.1",
xlab='Mortar Formulation',
ylab='Strength'
)
\noindent This plot also suggests that the modified mortar formulation has a lower mean strength than the unmodified formulation. Note that I assigned the return value from the boxplot function to a variable bp.returned. Several graphics functions return data used in the graphical construction. The boxplot function returns the summary statistics used to create the plots:
print(bp.returned)
Statistics and Distributions
The box plot was a nice way to illustrate the distribution of our data, but we could also use the summary function to get summary statistics:
summary(Table2.1)
Other statistics of interest, such as standard deviation (or variance) and empirical quantiles, are easy to collect. I'll define a new data.frame with some statistics:
dstat <- data.frame(
'Modified'=c(
mean(ModifiedMortar),
var(ModifiedMortar),
sd(ModifiedMortar),
length(ModifiedMortar)
),
'Unmodified'=c(
mean(UnmodifiedMortar),
var(UnmodifiedMortar),
sd(UnmodifiedMortar),
length(UnmodifiedMortar)
),
row.names=c('mean', 'variance', 'sd','obs')
)
\noindent The dstat object then contains the following:
kable(dstat) %>%
kable_styling(
bootstrap_options="striped",
full_width=FALSE
)
If we want to create a pooled estimator of variance $S_p^2$ we can do so manually:
n1 <- dstat['obs','Modified']
n2 <- dstat['obs','Unmodified']
S1sq <- dstat['var','Modified']
S2sq <- dstat['var','Unmodified']
pooled.var <- ((n1 - 1) * S1sq + (n2 - 1) * S2sq) / (n1 + n2 - 2)
print(pooled.var)
\noindent This calculation matches that of the text. Similarly, we can create the $t$-statistic to test the hypotheses $$
H_0: \mu_\text{Modified} = \mu_\text{Unmodified}
\quad\text{vs.}\quad
H_1: \mu_\text{Modified} \ne \mu_\text{Unmodified}~,
$$
ybar1 <- dstat['mean', 'Modified']
ybar2 <- dstat['mean', 'Unmodified']
t0 <- (ybar1 - ybar2) / sqrt(pooled.var * (1/n1 + 1/n2))
print(t0)
\noindent We can evaluate the statistic in two ways. First, we can use the critical region approach and find $t_{0.025, 18}$ using the qt function. Note that in R, all of the quantile functions give the area to the left (the lower area), but this text and many others take the symbol $t_{\alpha, n}$ to mean the critical point where the area to the right is $\alpha$. To adjust for this we give the argument lower.tail=FALSE.
critical.t <- qt(0.025, 18, lower.tail=FALSE)
print(critical.t)
\noindent Because $|t_0| > t_{0.025, 18}$ we would reject $H_0$. The other option is to find the $p$-value using the $pt$ function:
2 * (1 - pt(abs(t0), 18))
\noindent Since the $p$-value was less that $0.05$ we would again reject $H_0$.
That was a lot of typing to just perform a two-sample $t$-test! A shorter path would be to use R's t.test function. This requires us to reshape our data into a column for strength and a column that indicates if the observation was for the modified or unmodified formulation:
df <- data.frame(
'Strength'=c(ModifiedMortar, UnmodifiedMortar),
'Mortar'=factor(c(rep('Modified', 10), rep('Unmodified', 10)))
)
t.test(Strength ~ Mortar, data=df, paired=FALSE, var.equal=TRUE)
\noindent The t.test function is giving us a lot of information: the statistic, the degrees of freedom, the $p$-value, a confidence interval at 95%, and estimated means for the modified and unmodified mortars.
Section 2.4.3 includes the normal probability plots of tension bond strength in Figure 2.11. It's easiest to construct separate plots for each data set as follows:
op <- par(mfrow=c(1,2))
qqnorm(ModifiedMortar, main="Modified Mortar")
qqline(ModifiedMortar, col='orange')
qqnorm(UnmodifiedMortar, main="Unmodified Mortar")
qqline(UnmodifiedMortar, col='blue')
par(op)
\noindent Visual inspection doesn't suggest any strong violation of the assumption. Note that par is setting the graphical parameter mfrow so that we can plot two graphs next to each other. When setting with par the original value is returned and stored in op, which we then reset at the end.
detach(Table2.1) # Clean up since we're done with this data.
Power and Sample Size
For power and sample size calculations we'll use the pwr package [@pwr]. For the $t$-test we simply leave out the argument we want calculated. Recall that $$
d = \frac{|\mu_1 - \mu_2|}{\sigma}
$$
so for the mortar experiment we'd want $d=2$, then if each group has 10 observations we find the power:
library(pwr)
pwr.t.test(n=10, d=2, sig.level=0.05, type='two.sample', alternative='two.sided')
\noindent Here, our power was 98.8%. How many observations would we have needed to get the power of our test to 90%?
pwr.t.test(power=0.9, d=2, sig.level=0.05, type='two.sample', alternative='two.sided')
\noindent If we were planning the mortar experiment we'd want to make sure we had about 7 observations of each mortar to achieve approximately 90% power.
Unequal Variance --- Example 2.1
In Example 2.1 a study of the normalized fluorescence (after two hours) was undertaken for nerve tissue from 12 mice and muscle tissue from 12 mice. Note that our data is not paired. We wish to test $$
H_0: \mu_\text{Nerve} = \mu_\text{Muscle}
\quad\text{vs.}\quad
H_A: \mu_\text{Nerve} > \mu_\text{Muscle}~.
$$
\noindent The data from this experiment is presented in Table 2.3.
kable(Table2.3) %>% kable_styling()
\noindent Summary statistics for the data are presented below:
summary(Table2.3[,c('Nerve','Muscle')])
\noindent The summary statistics left out the standard deviations, which we can calculate below:
c('Nerve'=sd(Table2.3$Nerve), 'Muscle'=sd(Table2.3$Muscle))
\noindent Those estimates are strikingly different. Examining the normal probability plots we find a very distinct difference in slope, suggesting that the variances are unequal.
op <- par(mfrow=c(1,2))
qqnorm(Table2.3$Nerve, main="Nerve")
qqline(Table2.3$Nerve, col='orange')
qqnorm(Table2.3$Muscle, main="Muscle")
qqline(Table2.3$Muscle, col='blue')
par(op)
To test our hypothesis we again use the t.test procedure. To do the one-sided alternative test we use the argument alternative='greater'. Note that the order of the factor Tissue is Muscle then Nerve, so greater is testing that Nerve`'s mean is greater than Muscle's mean.
df <- data.frame(
'Flourescence'=c(Table2.3$Muscle, Table2.3$Nerve),
'Tissue'=factor(c(rep('Muscle', 12), rep('Nerve', 12)))
)
t.test(Flourescence ~ Tissue, data=df, paired=FALSE, var.equal=FALSE, alternative='greater')
Paired Comparisons
In Section 2.5.1 the hardness testing experiment is introduced. In this experiment, the depth of two different tips in samples of metal is measured. Note that each tip was tested on the same piece of metal, so our data is paired. The data for this experiment is given in Table 2.6.
kable(Table2.6) %>% kable_styling()
\noindent To perform the paired $t$-test we do the following:
df <- data.frame(
'Depth'=c(Table2.6$Tip1, Table2.6$Tip2),
'Tip'=factor(c(rep('Tip1', 10), rep('Tip2', 10)))
)
t.test(Depth ~ Tip, data=df, paired=TRUE)
\noindent With a $p$-value of 0.8 we certainly don't have evidence to reject our null hypothesis that the means are equal.
Notes
There are many other functions in R that handle a lot of the analysis we might want to do. A non-exhaustive list follows:
binom.test - Exact binomial test.chisq.test - Pearson's $\chi^2$-test for count data.krustkal.test - Kruskal-Wallis rank-sum test.prop.test - A one or two sample test of proportions.quantile - Get the empirical quantiles of a set of data.shapiro.test - Tests for normality of the data.summary - Produce a summary of data or other object.t.test - The two sample $t$-test.table - Tabulate data. Note that if you summary a table then a $\chi^2$ test of independence is performed.var.test - Conduct the two sample $F$-test for population variancewilcox.test - Test and construct confidence interval for the median.
References
ehassler/MontgomeryDAE documentation built on March 20, 2021, 7:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) library(kableExtra) library(MontgomeryDAE) source('_format.R')
Chapter 2 begins with several methods of visualization that are easily reproducible within R. These visualizations are based on the data in Table 2.1 which records the results of an experiment where a researcher was comparing the strength of a new mortar formulation to that of the original unmodified formulation.
kable(Table2.1) %>% kable_styling()
In R, the object Table2.1 is a data.frame that has a column ModifiedMortar and a column UnmodifiedMortar. The structure of this object is given below.
str(Table2.1)
\noindent We can use the R command attach to take the columns of the data.frame and make them accessible by name (e.g. we can use ModifiedMortar instead of Table2.1$ModifiedMortar):
attach(Table2.1) all(ModifiedMortar == Table2.1$ModifiedMortar)
\noindent When we're all done with the data we should call detach(Table2.1).
With the data loaded let's move on to some of the visualizations in Chapter 2.
Visualizations
We begin with the dot plot of Figure 2.1:
stripchart( list( 'Modified'=ModifiedMortar, 'Unmodified'=UnmodifiedMortar ), main="Dot Diagram for Table 2.1", xlab='Strength' )
\noindent The graph seems to suggest that the modified formulation has a lower average strength than that of the unmodified formulation.
Figure 2.2 in the text is a histogram of 200 observations of metal recovery for which no data is given. Though we lack the full data set, let's construct a histogram of the 20 data points from combining the modified and unmodified mortar formulation data:
data <- c(ModifiedMortar, UnmodifiedMortar) hist(data)
\noindent Though it was easy to construct, the histogram is not giving us a lot of information in this case due to the fact that we have such a small amount of data.
The dot plot of Figure 2.1 can become very difficult to interpret when the number of observations becomes large. In these cases we can use a box plot which is constructed based on summary statistics of the data. The box plots of Figure 2.3 are presented next.
bp.returned <- boxplot( Table2.1, main="Box plots for Table 2.1", xlab='Mortar Formulation', ylab='Strength' )
\noindent This plot also suggests that the modified mortar formulation has a lower mean strength than the unmodified formulation. Note that I assigned the return value from the boxplot function to a variable bp.returned. Several graphics functions return data used in the graphical construction. The boxplot function returns the summary statistics used to create the plots:
print(bp.returned)
Statistics and Distributions
The box plot was a nice way to illustrate the distribution of our data, but we could also use the summary function to get summary statistics:
summary(Table2.1)
Other statistics of interest, such as standard deviation (or variance) and empirical quantiles, are easy to collect. I'll define a new data.frame with some statistics:
dstat <- data.frame( 'Modified'=c( mean(ModifiedMortar), var(ModifiedMortar), sd(ModifiedMortar), length(ModifiedMortar) ), 'Unmodified'=c( mean(UnmodifiedMortar), var(UnmodifiedMortar), sd(UnmodifiedMortar), length(UnmodifiedMortar) ), row.names=c('mean', 'variance', 'sd','obs') )
\noindent The dstat object then contains the following:
kable(dstat) %>% kable_styling( bootstrap_options="striped", full_width=FALSE )
If we want to create a pooled estimator of variance $S_p^2$ we can do so manually:
n1 <- dstat['obs','Modified'] n2 <- dstat['obs','Unmodified'] S1sq <- dstat['var','Modified'] S2sq <- dstat['var','Unmodified'] pooled.var <- ((n1 - 1) * S1sq + (n2 - 1) * S2sq) / (n1 + n2 - 2) print(pooled.var)
\noindent This calculation matches that of the text. Similarly, we can create the $t$-statistic to test the hypotheses $$ H_0: \mu_\text{Modified} = \mu_\text{Unmodified} \quad\text{vs.}\quad H_1: \mu_\text{Modified} \ne \mu_\text{Unmodified}~, $$
ybar1 <- dstat['mean', 'Modified'] ybar2 <- dstat['mean', 'Unmodified'] t0 <- (ybar1 - ybar2) / sqrt(pooled.var * (1/n1 + 1/n2)) print(t0)
\noindent We can evaluate the statistic in two ways. First, we can use the critical region approach and find $t_{0.025, 18}$ using the qt function. Note that in R, all of the quantile functions give the area to the left (the lower area), but this text and many others take the symbol $t_{\alpha, n}$ to mean the critical point where the area to the right is $\alpha$. To adjust for this we give the argument lower.tail=FALSE.
critical.t <- qt(0.025, 18, lower.tail=FALSE) print(critical.t)
\noindent Because $|t_0| > t_{0.025, 18}$ we would reject $H_0$. The other option is to find the $p$-value using the $pt$ function:
2 * (1 - pt(abs(t0), 18))
\noindent Since the $p$-value was less that $0.05$ we would again reject $H_0$.
That was a lot of typing to just perform a two-sample $t$-test! A shorter path would be to use R's t.test function. This requires us to reshape our data into a column for strength and a column that indicates if the observation was for the modified or unmodified formulation:
df <- data.frame( 'Strength'=c(ModifiedMortar, UnmodifiedMortar), 'Mortar'=factor(c(rep('Modified', 10), rep('Unmodified', 10))) ) t.test(Strength ~ Mortar, data=df, paired=FALSE, var.equal=TRUE)
\noindent The t.test function is giving us a lot of information: the statistic, the degrees of freedom, the $p$-value, a confidence interval at 95%, and estimated means for the modified and unmodified mortars.
Section 2.4.3 includes the normal probability plots of tension bond strength in Figure 2.11. It's easiest to construct separate plots for each data set as follows:
op <- par(mfrow=c(1,2)) qqnorm(ModifiedMortar, main="Modified Mortar") qqline(ModifiedMortar, col='orange') qqnorm(UnmodifiedMortar, main="Unmodified Mortar") qqline(UnmodifiedMortar, col='blue') par(op)
\noindent Visual inspection doesn't suggest any strong violation of the assumption. Note that par is setting the graphical parameter mfrow so that we can plot two graphs next to each other. When setting with par the original value is returned and stored in op, which we then reset at the end.
detach(Table2.1) # Clean up since we're done with this data.
Power and Sample Size
For power and sample size calculations we'll use the pwr package [@pwr]. For the $t$-test we simply leave out the argument we want calculated. Recall that $$
d = \frac{|\mu_1 - \mu_2|}{\sigma}
$$
so for the mortar experiment we'd want $d=2$, then if each group has 10 observations we find the power:
library(pwr) pwr.t.test(n=10, d=2, sig.level=0.05, type='two.sample', alternative='two.sided')
\noindent Here, our power was 98.8%. How many observations would we have needed to get the power of our test to 90%?
pwr.t.test(power=0.9, d=2, sig.level=0.05, type='two.sample', alternative='two.sided')
\noindent If we were planning the mortar experiment we'd want to make sure we had about 7 observations of each mortar to achieve approximately 90% power.
Unequal Variance --- Example 2.1
In Example 2.1 a study of the normalized fluorescence (after two hours) was undertaken for nerve tissue from 12 mice and muscle tissue from 12 mice. Note that our data is not paired. We wish to test $$ H_0: \mu_\text{Nerve} = \mu_\text{Muscle} \quad\text{vs.}\quad H_A: \mu_\text{Nerve} > \mu_\text{Muscle}~. $$
\noindent The data from this experiment is presented in Table 2.3.
kable(Table2.3) %>% kable_styling()
\noindent Summary statistics for the data are presented below:
summary(Table2.3[,c('Nerve','Muscle')])
\noindent The summary statistics left out the standard deviations, which we can calculate below:
c('Nerve'=sd(Table2.3$Nerve), 'Muscle'=sd(Table2.3$Muscle))
\noindent Those estimates are strikingly different. Examining the normal probability plots we find a very distinct difference in slope, suggesting that the variances are unequal.
op <- par(mfrow=c(1,2)) qqnorm(Table2.3$Nerve, main="Nerve") qqline(Table2.3$Nerve, col='orange') qqnorm(Table2.3$Muscle, main="Muscle") qqline(Table2.3$Muscle, col='blue') par(op)
To test our hypothesis we again use the t.test procedure. To do the one-sided alternative test we use the argument alternative='greater'. Note that the order of the factor Tissue is Muscle then Nerve, so greater is testing that Nerve`'s mean is greater than Muscle's mean.
df <- data.frame( 'Flourescence'=c(Table2.3$Muscle, Table2.3$Nerve), 'Tissue'=factor(c(rep('Muscle', 12), rep('Nerve', 12))) ) t.test(Flourescence ~ Tissue, data=df, paired=FALSE, var.equal=FALSE, alternative='greater')
Paired Comparisons
In Section 2.5.1 the hardness testing experiment is introduced. In this experiment, the depth of two different tips in samples of metal is measured. Note that each tip was tested on the same piece of metal, so our data is paired. The data for this experiment is given in Table 2.6.
kable(Table2.6) %>% kable_styling()
\noindent To perform the paired $t$-test we do the following:
df <- data.frame( 'Depth'=c(Table2.6$Tip1, Table2.6$Tip2), 'Tip'=factor(c(rep('Tip1', 10), rep('Tip2', 10))) ) t.test(Depth ~ Tip, data=df, paired=TRUE)
\noindent With a $p$-value of 0.8 we certainly don't have evidence to reject our null hypothesis that the means are equal.
Notes
There are many other functions in R that handle a lot of the analysis we might want to do. A non-exhaustive list follows:
binom.test- Exact binomial test.chisq.test- Pearson's $\chi^2$-test for count data.krustkal.test- Kruskal-Wallis rank-sum test.prop.test- A one or two sample test of proportions.quantile- Get the empirical quantiles of a set of data.shapiro.test- Tests for normality of the data.summary- Produce a summary of data or other object.t.test- The two sample $t$-test.table- Tabulate data. Note that if yousummarya table then a $\chi^2$ test of independence is performed.var.test- Conduct the two sample $F$-test for population variancewilcox.test- Test and construct confidence interval for the median.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.