Before we can actually use the package, we’ll need to load it. We’ll
also load the dplyr package, which we’ll talk about
shortly. There are also a few other packages that we’ll load. These
could all be loaded by loading the tidyverse package
instead.
Calculating Single-Point Basis Value
The simplest thing that you will likely do is to calculate a basis
value based of a set of numbers that you consider as unstructured data.
An example of this would be calculating the B-Basis of the
RTD warp tension (WT) data.
There are a number of diagnostic tests that we should run before
actually calculating a B-Basis value. We’ll talk about those later, but
for now, let’s just get right to checking how the data are distributed
and calculating the B-Basis.
We’ll use an Anderson–Darling test to check if the data are normally
distributed. The cmstatr package provides the function
anderson_darling_normal and related functions for other
distributions. We can run an Anderson–Darling test for normality on the
warp tension RTD data as follows. We’ll perform this test on the
normalized strength.
norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
anderson_darling_normal(strength.norm)
#>
#> Call:
#> anderson_darling_normal(data = ., x = strength.norm)
#>
#> Distribution: Normal ( n = 28 )
#> Test statistic: A = 0.3805995
#> OSL (p-value): 0.3132051 (assuming unknown parameters)
#> Conclusion: Sample is drawn from a Normal distribution ( alpha = 0.05 )
Now that we know that this data follows a normal distribution (since
the observed significance level (OSL) of the Anderson–Darling test is
greater than \(0.05\)), we can proceed
to calculate a basis value based based on the assumption of normally
distributed data. The cmstatr package provides the function
basis_normal as well as related functions for other
distributions. By default, the B-Basis value is calculated, but other
population proportions and confidence bounds can be specified (for
example, specify p = 0.99, conf = 0.99 for A-Basis).
norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
basis_normal(strength.norm)
#> `outliers_within_batch` not run because parameter `batch` not specified
#> `between_batch_variability` not run because parameter `batch` not specified
#>
#> Call:
#> basis_normal(data = ., x = strength.norm)
#>
#> Distribution: Normal ( n = 28 )
#> B-Basis: ( p = 0.9 , conf = 0.95 )
#> 129.9583
We see that the calculated B-Basis is \(129.96\). We also see two messages issued
by the cmstatr package. These messages relate to the
automated diagnostic tests performed by the basis calculation functions.
In this case we see messages that two of the diagnostic tests were not
performed because we didn’t specify the batch of each observation. The
batch is not required for calculating single-point basis values, but it
is required for performing batch-to-batch variability and within-batch
outlier diagnostic tests.
The basis_normal function performs the following
diagnostic tests by default:
- Within batch outliers using
maximum_normed_residual()
- Between batch variability using
ad_ksample()
- Outliers using
maximum_normed_residual()
- Normality of data using
anderson_darling_normal()
There are two ways that we can deal with the two messages that we
see. We can pass in a column that specifies the batch for each
observation, or we can override those two diagnostic tests so that
cmstatr doesn’t run them.
To override the two diagnostic tests, we set the argument
override to a list of the names of the diagnostic tests
that we want to skip. The names of the diagnostic tests that were not
run are shown between back-ticks (`) in the message. Our call to
basis_normal() would be updated as follows:
norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
basis_normal(strength.norm,
override = c("outliers_within_batch",
"between_batch_variability"))
#>
#> Call:
#> basis_normal(data = ., x = strength.norm, override = c("outliers_within_batch",
#> "between_batch_variability"))
#>
#> Distribution: Normal ( n = 28 )
#> The following diagnostic tests were overridden:
#> `outliers_within_batch`,
#> `between_batch_variability`
#> B-Basis: ( p = 0.9 , conf = 0.95 )
#> 129.9583
Obviously, you should be cautious about overriding the diagnostic
tests. There are certainly times when it is appropriate to do so, but
sound engineering judgment is required.
The better approach would be to specify the batch. This can be done
as follows. We’ll store the result in the variable
b_basis_wt_rtd for reasons that will become clear
later.
b_basis_wt_rtd <- norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
basis_normal(strength.norm, batch)
#> Warning: `between_batch_variability` failed: Anderson-Darling k-Sample test
#> indicates that batches are drawn from different distributions
Now that batch is specified, we see that one of the diagnostic tests
actually fails: the Anderson–Darling k-Sample test shows that the
batches are not drawn from the same (unspecified) distribution. We can
interrogate the failing test by accessing the
diagnostic_obj element of the return value from
basis_normal(). This contains elements for each of the
diagnostic tests. We can access the
between_batch_variability result as follows:
b_basis_wt_rtd$diagnostic_obj$between_batch_variability
#>
#> Call:
#> ad_ksample(x = x, groups = batch, alpha = 0.025)
#>
#> N = 28 k = 3
#> ADK = 6.65 p-value = 0.0025892
#> Conclusion: Samples do not come from the same distribution (alpha = 0.025 )
We could have also run the failing diagnostic test directly as
follows:
norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
ad_ksample(strength.norm, batch)
#>
#> Call:
#> ad_ksample(data = ., x = strength.norm, groups = batch)
#>
#> N = 28 k = 3
#> ADK = 6.65 p-value = 0.0025892
#> Conclusion: Samples do not come from the same distribution (alpha = 0.025 )
For the Anderson–Darling k-Sample test, \(\alpha=0.025\) is normally used. In this
case the p-value is \(p=0.0026\), so it
is no where near \(\alpha\) (note the
number of decimal places).
We can plot the distribution of this data and make a judgment call
about whether to continue.
norm_data %>%
filter(test == "WT" & condition == "RTD") %>%
group_by(batch) %>%
ggplot(aes(x = strength.norm, color = batch)) +
stat_normal_surv_func() +
stat_esf() +
ggtitle("Distribution of Data For Each Batch")
We can also run the other diagnostic test by themselves. These are
described in more detail in the following sections.
Batch-to-Batch Distribution
Next, we will use the Anderson–Darling k-Sample test to check that
each batch comes from the same distribution within each condition. We
can use the ad_ksample function from cmstatr
to do so. Once again, we’ll use nest/unnest
and glance to do so.
norm_data %>%
filter(test == "FC") %>%
group_by(condition) %>%
nest() %>%
mutate(adk = map(data, ~ad_ksample(data = .x,
x = strength.norm,
groups = batch)),
tidied = map(adk, glance)) %>%
select(-c(data, adk)) %>% # remove unneeded columns
unnest(tidied)
#> # A tibble: 5 × 8
#> # Groups: condition [5]
#> condition alpha n k sigma ad p reject_same_dist
#> <chr> <dbl> <int> <int> <dbl> <dbl> <dbl> <lgl>
#> 1 CTD 0.025 18 3 0.944 1.76 0.505 FALSE
#> 2 RTD 0.025 18 3 0.944 1.03 0.918 FALSE
#> 3 ETD 0.025 18 3 0.944 0.683 0.997 FALSE
#> 4 ETW 0.025 18 3 0.944 0.93 0.954 FALSE
#> 5 ETW2 0.025 19 3 0.951 1.74 0.513 FALSE
For all conditions, the Anderson–Darling k-Sample test fails to
reject the hypothesis that each batch comes from the same (unspecified)
distribution. We can thus proceed to pooling the data.
Checking for Outliers Within Each Condition
Just as we did when checking for outlier within each condition and
each batch, we can pool all the batches (within each condition) and
check for outliers within each condition.
norm_data %>%
filter(test == "FC") %>%
group_by(condition) %>%
nest() %>%
mutate(mnr = map(data, ~maximum_normed_residual(data = .x,
x = strength.norm)),
tidied = map(mnr, glance)) %>%
select(-c(mnr, data)) %>% # remove unneeded columns
unnest(tidied)
#> # A tibble: 5 × 5
#> # Groups: condition [5]
#> condition mnr alpha crit n_outliers
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 CTD 2.38 0.05 2.65 0
#> 2 RTD 2.06 0.05 2.65 0
#> 3 ETD 2.05 0.05 2.65 0
#> 4 ETW 2.34 0.05 2.65 0
#> 5 ETW2 2.07 0.05 2.68 0
We find no outliers, so we can continue.
Investigation Conditions
When multiple conditions were tested, it’s usually useful to view
some basic summary statistics for each condition before proceeding. The
condition_summary function can be used for this. You can
pass a data.frame with the data and the name of the
condition variable to generate such summary statistics.
norm_data %>%
filter(test == "FC") %>%
condition_summary(strength.norm, condition, "RTD")
#> condition n mean mean_fraction
#> 1 CTD 18 96.41519 1.0877571
#> 2 RTD 18 88.63668 1.0000000
#> 3 ETD 18 75.40668 0.8507390
#> 4 ETW 18 58.28060 0.6575223
#> 5 ETW2 19 51.87082 0.5852071
Pooling Across Environments
Often it is desirable to pool data across several environments. There
are two methods for doing so: “pooled standard deviation” and “pooled
CV” (CV is an abbreviation for Coefficient of Variation).
First, we will check for equality of variance among the conditions.
We will do so using Levene’s test. The cmstatr package
provides the function levene_test to do so.
norm_data %>%
filter(test == "FC") %>%
levene_test(strength.norm, condition)
#>
#> Call:
#> levene_test(data = ., x = strength.norm, groups = condition)
#>
#> n = 91 k = 5
#> F = 5.260731 p-value = 0.0007727083
#> Conclusion: Samples have unequal variance ( alpha = 0.05 )
The result from Levene’s test indicates that the variance for each
condition is not equal. This indicates that the data cannot be pooled
using the “pooled standard deviation” method.
We can check if the data can be pooled using the “pooled CV” method.
We’ll start by normalizing the data from each group to the group’s mean.
The cmstatr package provides the function
normalize_group_mean for this purpose.
norm_data %>%
filter(test == "FC") %>%
mutate(
strength_norm_group = normalize_group_mean(strength.norm, condition)) %>%
levene_test(strength_norm_group, condition)
#>
#> Call:
#> levene_test(data = ., x = strength_norm_group, groups = condition)
#>
#> n = 91 k = 5
#> F = 1.839645 p-value = 0.1285863
#> Conclusion: Samples have equal variances ( alpha = 0.05 )
The Levene’s test thus shows the variances of the pooled data are
equal. We can move on to performing an Anderson–Darling test for
normality on the pooled data.
norm_data %>%
filter(test == "FC") %>%
mutate(
strength_norm_group = normalize_group_mean(strength.norm, condition)) %>%
anderson_darling_normal(strength_norm_group)
#>
#> Call:
#> anderson_darling_normal(data = ., x = strength_norm_group)
#>
#> Distribution: Normal ( n = 91 )
#> Test statistic: A = 0.3619689
#> OSL (p-value): 0.3812268 (assuming unknown parameters)
#> Conclusion: Sample is drawn from a Normal distribution ( alpha = 0.05 )
The Anderson–Darling test indicates that the pooled data is drawn
from a normal distribution, so we can continue with calculating basis
values using the “pooled CV” method.
norm_data %>%
filter(test == "FC") %>%
basis_pooled_cv(strength.norm, condition, batch)
#>
#> Call:
#> basis_pooled_cv(data = ., x = strength.norm, groups = condition,
#> batch = batch)
#>
#> Distribution: Normal - Pooled CV ( n = 91, r = 5 )
#> B-Basis: ( p = 0.9 , conf = 0.95 )
#> CTD 85.09241
#> ETD 66.55109
#> ETW 51.43626
#> ETW2 45.81318
#> RTD 78.2274
The conditions listed in the output above are in alphabetical order.
This probably isn’t what you want. Instead, you probably want the
conditions listed in a certain order. This can be done by ordering the
data first as demonstrated below. You’re probably just do this one in at
the start of your analysis. In the example below, we’ll store the result
in the variable basis_res before printing it.
basis_res <- norm_data %>%
mutate(condition = ordered(condition,
c("CTD", "RTD", "ETD", "ETW", "ETW2"))) %>%
filter(test == "FC") %>%
basis_pooled_cv(strength.norm, condition, batch)
basis_res
#>
#> Call:
#> basis_pooled_cv(data = ., x = strength.norm, groups = condition,
#> batch = batch)
#>
#> Distribution: Normal - Pooled CV ( n = 91, r = 5 )
#> B-Basis: ( p = 0.9 , conf = 0.95 )
#> CTD 85.09241
#> RTD 78.2274
#> ETD 66.55109
#> ETW 51.43626
#> ETW2 45.81318
The summary statistics that we computed for each condition earlier
can also be generated using the basis object returned by
basis_pooled_cv and related functions.
basis_res %>%
condition_summary("RTD")
#> condition n mean mean_fraction basis basis_fraction
#> 1 CTD 18 96.41519 1.0877571 85.09241 1.0877571
#> 2 RTD 18 88.63668 1.0000000 78.22740 1.0000000
#> 3 ETD 18 75.40668 0.8507390 66.55109 0.8507390
#> 4 ETW 18 58.28060 0.6575223 51.43626 0.6575223
#> 5 ETW2 19 51.87082 0.5852071 45.81318 0.5856412
Equivalency
Eventually, once you’ve finished calculating all your basis values,
you’ll probably want to set specification requirements or evaluate
site/process equivalency. cmstatr has functionality to do
both.
Let’s say that you want to develop specification limits for fill
compression that you’re going to put in your material specification. You
can do this as follows:
carbon.fabric.2 %>%
filter(test == "FC" & condition == "RTD") %>%
equiv_mean_extremum(strength, n_sample = 5, alpha = 0.01)
#>
#> Call:
#> equiv_mean_extremum(df_qual = ., data_qual = strength, n_sample = 5,
#> alpha = 0.01)
#>
#> For alpha = 0.01 and n = 5
#> ( k1 = 3.071482 and k2 = 1.142506 )
#> Min Individual Sample Mean
#> Thresholds: 69.89842 82.16867
If you’re determining equivalency limits for modulus, a different
approach is generally used so that bilateral limits are set.
cmstatr can do this as well, using the function
equiv_change_mean.