## Comparing Two Sample Means in R

One can easily compare two sample means in R, as in R language all the classical tests are available in the package stats. There are different comparison tests such as (i) one sample mean test, (ii) two independent sample means test, and (iii) dependent sample test. When population standard deviation is known, or sample size (number of observations in the sample) is large enough ($n\ge 30), tests related to normal distribution are performed. ### Data for Two Sample Means Consider the following data set on the “latent heat of the fusion of ice (cal/gm)” from Rice, 1995. Let us draw boxplots to make a comparison between two these two methods. The comparison will help in checking the assumption of the independent two-sample test. Note that one can read the data using the scan() function, create vectors, or even read the above data from data files such as *.txt and *.csv. In this tutorial, we assume vectors$A$and$B$for method A and method B. A = c(79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97, 80.05, 80.03, 80.02, 80.00, 80.02) B = c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97) ### Draw a Boxplot of Samples Let us draw boxplots for each method that indicate the first group tends to give higher results than the second one. boxplot(A, B) ### Comparing Two Sample Means in R using t.test() Function The unpaired t-test (independent two-sample test) for the equality of the means can be done using the function t.test() in R Language. t.test(A, B) From the results above, one can see that the p-value = 0.006939 is less than 0.05 (level of significance) which means that on average both methods are statistically different from each other with reference to latent heat of fusion of ice. ### Testing the Equality of Variances of Samples Note that, the R language does not assume the equality of variances in the two samples. However, the F-test can be used to check/test the equality in the variances, provided that the two samples are from normal populations. var.test(A, B) From the above results, there is no evidence that the variances of both samples are statistically significant, as the p-value is greater than the 0.05 level of significance. It means that one can use the classical t-test that assumes the equality of the variances. t.test(A, B, var.equa. = TRUE) ## Output Welch Two Sample t-test data: A and B t = 3.2499, df = 12.027, p-value = 0.006939 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.01385526 0.07018320 sample estimates: mean of x mean of y 80.02077 79.97875  https://rfaqs.com https://gmstat.com ## Statistical Power Analysis in R: A Comprehensive Guide ### Introduction to Power Analysis The post is about statistical power analysis in R. First, define the meaning of power in statistics. The power is the probability ($1-\beta$) of detecting an effect given that the effect is here. Power is the probability of correctly rejecting the null hypothesis when it is false. Suppose, a simple study of a drug-A and a placebo. Let the drug be truly effective. The power is the probability of finding a difference between two groups (drug-A and placebo group). Imagine that a power of$1-\beta=0.8$(having a power of 0.8 means that 80% of the time, there will be statistically significant differences between the drug-A and the placebo group, whereas there are 20% of the time, the statistically significant effect will not be obtained between two groups). Also, note that this study was conducted many times. Therefore, the probability of a Type-II error is$\beta=0.2$. ### One-Sample Power The following plot is for a one-sample one-tailed greater than t-test. In the graph below, let the null hypothesis$H_0:\mu = \mu_0$be true, and the test statistic$t$follows the null distribution indicated by the hashed area. Under the specific alternative hypothesis,$H_1:\mu = \mu_1$, the test statistic$t$follows the distribution shown by solid area. The$\alpha$is the probability of making a type-I error (that is rejecting$H_0$when it is true), and the “crit. Val” is the location of the$t_{crit}$value associated with$H_0$on the scale of the data. The rejection region is the area under$H_0$at least as far as$crit. val.” is from $\mu_0$.

The test’s power ($1-\beta$) is the green area, the area under $H_1$ in the rejection region. A type-II error is made when $H_1$ is true, but we fail to reject $H_0$ in the red region.

### Type-II Error and Power Analysis in R

#One Sample Power

x <- seq(-4, 4, length = 1000)
hx <- dnorm(x, mean = 0, sd = 1)

plot(x, hx, type = "n", xlim = c(-4, 8), ylim = c(0, 0.5),
main = expression (paste("Type-II Error (", beta, ") and Power (", 1 - beta, ")")),
axes = FALSE)

# one-tailed shift
shift = qnorm (1 - 0.05, mean=0, sd = 1 )*1.7
xfit2 = x + shift
yfit2 = dnorm(xfit2, mean=shift, sd = 1 )

axis (1, at = c(-qnorm(0.05), 0, shift), labels = expression("crit. val.", mu[0], mu[1]))
axis(1, at = c(-4, 4 + shift), labels = expression(-infinity, infinity),
lwd = 1, lwd.tick = FALSE)

# The alternative hypothesis area
# the red - underpowered area

lb <- min(xfit2)               # lower bound
ub <- round(qnorm(0.95), 2)    # upper bound
col1 = "#CC2222"

i <- xfit2 >= lb & xfit2 <= ub
polygon(c(lb, xfit2[i], ub), c(0, yfit2[i],0), col = col1)

# The green area where the power is
col2 = "#22CC22"
i <- xfit2 >= ub
polygon(c(ub, xfit2[i], max(xfit2)), c(0, yfit2[i], 0), col = col2)

# Outline the alternative hypothesis
lines(xfit2, yfit2, lwd = 2)

# Print null hypothesis area
col_null = "#AAAAAA"
polygon (c(min(x), x, max(x)), c(0, hx, 0), col = col_null,
lwd = 2, density = c(10, 40), angle = -45, border = 0)

lines(x, hx, lwd = 2, lty = "dashed", col=col_null)

axis(1, at = (c(ub, max(xfit2))), labels = c("", expression(infinity)), col = col2,
lwd = 1, lwd.tick = FALSE)

#Legend
legend("topright", inset = 0.015, title = "Color",
c("Null Hypothesis", "Type-II error", "Power"), fill = c(col_null, col1, col2),
angle = -45, density = c(20, 1000, 1000), horiz = FALSE)

abline(v=ub, lwd=2, col="#000088", lty = "dashed")
arrows(ub, 0.45, ub+1, 0.45, lwd=3, col="#008800")
arrows(ub, 0.45, ub-1, 0.45, lwd=3, col="#880000")

Online Quiz Website

Statistics and Data Analysis

## Mean Comparison Tests: Hypothesis Testing (One Sample and Two Sample)

Here we learn some basics about how to perform Mean Comparison Tests: hypothesis testing for one sample test, two-sample independent test, and dependent sample test. We will also learn how to find the p-values for a certain distribution such as t-distribution, and critical region values. We will also see how to perform one-tailed and two-tailed hypothesis tests.

### How to Perform One-Sample t-Test in R

A recent article in The Wall Street Journal reported that the 30-year mortgage rate is now less than 6%. A sample of eight small banks in the Midwest revealed the following 30-year rates (in percent)

At the 0.01 significance level (probability of type-I error), can we conclude that the 30-year mortgage rate for small banks is less than 6%?

#### Manual Calculations for One-Sample t-Test and Confidence Interval

One sample mean comparison test can be performed manually.

# Manual way
X <- c(4.8, 5.3, 6.5, 4.8, 6.1, 5.8, 6.2, 5.6)
xbar <- mean(X)
s <- sd(X)
mu = 6
n = length(X)
df = n - 1
tcal = (xbar - mu)/(s/sqrt(n) )
tcal
c(xbar - qt(0.995, df = df) * s/sqrt(n), xbar + qt(0.995, df = df) * s/sqrt(n))

#### Critical Values from t-Table

# Critical Value for Left Tail
qt(0.01, df = df, lower.tail = T)
# Critical Value for Right Tail
qt(0.99, df = df, lower.tail = T)
# Critical Vale for Both Tails
qt(0.995, df = df)

#### Finding p-Values

# p-value (altenative is less)
pt(tcal, df = df)
# p-value (altenative is greater)
1 - pt(tcal, df = df)
# p-value (alternative two tailed or not equal to)
2 * pt(tcal, df = df)

#### Performing One-Sample Confidence Interval and t-test Using Built-in Function

One can perform one sample mean comparison test using built-in functions available in the R Language.

# Left Tail test
t.test(x = X, mu = 6, alternative = c("less"), conf.level = 0.99)
# Right Tail test
t.test(x = X, mu = 6, alternative = c("greater"), conf.level = 0.99)
# Two Tail test
t.test(x = X, mu = 6, alternative = c("two.sided"), conf.level = 0.99)

### How to Perform two-Sample t-Test in R

Consider we have two samples stored in two vectors $X$ and $Y$ as shown in R code. We are interested in the Mean Comparison Test among two groups of people regarding (say) their wages in a certain week.

X = c(70, 82, 78, 70, 74, 82, 90)
Y = c(60, 80, 91, 89, 77, 69, 88, 82)

#### Manual Calculations for Two-Sample t-Test and Confidence Interval

The manual calculation for two sample t-tests as mean comparison test is as follows.

nx = length(X)
ny = length(Y)
xbar = mean(X)
sx = sd(X)
ybar = mean(Y)
sy = sd(Y)
df = nx + ny - 2
# Pooled Standard Deviation/ Variance
SP = sqrt( ( (nx-1) * sx^2 + (ny-1) * sy^2) / df )
tcal = (( xbar - ybar ) - 0) / (SP *sqrt(1/nx + 1/ny))
tcal
# Confidence Interval
LL <- (xbar - ybar) - qt(0.975, df)* sqrt((SP^2 *(1/nx + 1/ny) ))
UL <- (xbar - ybar) + qt(0.975, df)* sqrt((SP^2 *(1/nx + 1/ny) ))
c(LL, UL)

#### Finding p-values

# The p-value at the left-hand side of Critical Region
pt(tcal, df )
# The p-value for two-tailed Critical Region
2 * pt(tcal, df )
# The p-value at the right-hand side of Critical Region
1 - pt(tcal, df)

#### Finding Critical Values from t-Table

# Left Tail
qt(0.025, df = df, lower.tail = T)
# Right Tail
qt(0.975, df = df, lower.tail = T)
# Both tails
qt(0.05, df = df)

#### Performing Two-Sample Confidence Interval and T-test using Built-in Function

One can perform two sample mean comparison test using built-in functions in R Language.

# Left Tail test
t.test(X, Y, alternative = c("less"), var.equal = T)
# Right Tail test
t.test(X, Y, alternative = c("greater"), var.equal = T)
# Two Tail test
t.test(X, Y, alternative = c("two.sided"), var.equal = T)

Note if $X$ and $Y$ variables are from a data frame then perform the two-sample t-test using the formula symbol (~). Let’s first make the data frame from vectors $X$ and $$Y.

data <- data.frame(values = c(X, Y), group = c(rep("A", nx), rep("B", ny)))
t.test(values ~ group, data = data, alternative = "less", var.equal = T)
t.test(values ~ group, data = data, alternative = "greater", var.equal = T)
t.test(values ~ group, data = data, alternative = "two.side", var.equal = T)

To understand probability distributions functions in R click the link: Probability Distributions in R

MCQs in Statistics

## One-Way ANOVA in R: A Comprehensive Quide

In this post, we will learn about one-way ANOVA in R Language.

The two-sample t or z-test is used to compare two groups from the independent population. However, if there are more than two groups, One-Way ANOVA (analysis of variance) or its further versions can be used in R.

### Introduction to One-Way ANOVA

The statistical test statistic associated with ANOVA is the F-test (also called F-ratio). In the Anova procedure, an observed F-value is computed and then compared with a critical F-value derived from the relevant F-distribution. The F-value comes from a family of F-distribution defined by two numbers (the degrees of freedom). Note that the F-distribution cannot be negative as it is the ratio of variance and variances are always positive numbers.

The One-Way ANOVA is also known as one-factor ANOVA. It is the extension of the independent two-sample test for comparing means when there are more than two groups. The data in One-Way ANOVA is organized into several groups based on grouping variables (called factor variables too).

To compute the F-value, the ratio of “the variance between groups”,  and the “variance within groups” needs to be computed. The assumptions of ANOVA should also be checked before performing the ANOVA test. We will learn how to perform One-Way ANOVA in R.

### One-Way ANOVA in R

Suppose we are interested in finding the difference of miles per gallon based on number of the cylinders in an automobile; from the dataset “mtcars”

Let us get some basic insight into the data before performing the ANOVA.

# load and attach the data mtcars
attach(mtcars)
# see the variable names and initial observations
head(mtcars)

Let us find the means of each number of the cylinder group

print(model.tables(res, "means"), digits = 4)

Let us draw the boxplot of each group

boxplot(mpg ~ cyl, main="Boxplot", xlab="Number of Cylinders", ylab="mpg")

Now, to perform One-Way ANOVA in R using the aov( ) function. For example,

aov(mpg ~ cyl)

The variable “mpg” is continuous and the variable “cyl” is the grouping variable. From the output note the degrees of freedom under the variable “cyl”. It will be one. It means the results are not correct as the degrees of freedom should be two as there are three groups on “cyl”. In the mode (data type) of grouping variable required for ANOVA  should be the factor variable. For this purpose, the “cyl” variable can be converted to factor as

cyl <- as.factor(cyl)

Now re-issue the aov( ) function as

aov(mpg ~ cyl)

Now the results will be as required. To get the ANOVA table, use the summary( ) function as

summary(aov (mpg ~ cyl))

Let’s store the ANOVA results obtained from aov( ) in object say res

res <- aov(mpg ~ cyl)
summary(res)

### Post-Hoc Analysis (Multiple Pairwise Comparison)

Post-hoc tests or multiple-pairwise comparison tests help in finding out which groups differ (significantly) from one other and which do not. The post-hoc tests allow for multiple-pairwise comparisons without inflating the type-I error. To understand it, suppose the level of significance (type-I error) is 5%. Then the probability of making at least one Type-I error (assuming independence of three events), the maximum family-wise error rate will be

$1-(0.95 \times 0.95 \times 0.95) = 14.2%$

It will give the probability of having at least one FALSE alarm (type-I error).

To perform Tykey’s post-hoc test and plot the group’s differences in means from Tukey’s test.

# Tukey Honestly Significant Differences
TukeyHSD(res)
plot(TukeyHSD(res))

### Diagnostic Plots (Checking Model Assumptions)

The diagnostic plots can be used to check the assumption of heteroscedasticity, normality, and influential observations.

layout(matrix(c(1,2,3,4), 2,2))
plot(res)

#### Levene’s Test

To check the assumption of ANOVA, Levene’s test can be used. For this purpose leveneTest( ) function can be used which is available in the car package.

library(car)
leveneTest(res)

https://itfeature.com

https://gmstat.com