If \(X\) is normaly distributed with mean \(\mu\) and variance \(\sigma^2\), then the probability density function for a given draw from that distribution \(x\) is given by
\[ P(x) = \frac{1}{{\sigma \sqrt {2\pi } }}e^{{{ - \left( {x - \mu } \right)^2 } {2\sigma ^2 }}} \]
Further, the cumulative density function can be found by integration of this formula \[ P(X \leq x) = \frac{1}{{\sigma \sqrt {2\pi } }} \int_{-\inf}^x e^{{{ - \left( {x - \mu } \right)^2 } {2\sigma ^2 }}} \, dx \]
The short notation for this distribution is
\[ X \sim \mathcal{N}(\mu,\sigma^2) \]
More often than not, the distribution is unknown. That is the values of \(\mu\) and \(\sigma\) are unknown and must then be estimated from data.
\(\mu\) characterizes the center of the distribution, and is naturally estimated by the mean-value of the data-points (\(\bar{X}\)). \(\sigma\) reflects the spread around the mean, and is in a similar fashion estimated by the standard deviation (\(\hat{\sigma}\) or \(s\)).
\[ \hat{\mu} = \bar{X} = \sum_{i = 1}^n X_i / n = (X_1 + X_2 + X_3 + \ldots + X_n) / n \]
\[ \hat{\sigma}^2 = s^2 = \sum_{i = 1}^n (X_i - \bar{X})^2 / (n-1) = ((X_1 - \bar{X})^2 + (X_2 - \bar{X})^2 + (X_3 - \bar{X})^2 + \ldots + (X_n - \bar{X})^2 ) / (n-1) \]
Caffeine is a central neurvous system stimulant, which can have several positive- and negative effects. In this study, the level of activity is up for examination, and for this purpose a model-system; the runing activity of mice in a wheel reflected as the number of rounds pr minutes recorded over 7 minutes. Here a total of 249 mice from two species; one breed to be high in running performance, and one control, where given either Water, Gatorade or Red Bull followed by measuring their voluentary wheel run activiy. Of interest is the average actvity within each mouse type and for each caffeine type and how large the spread is.
The data is listed in a table with 249 rows (here, the first five samples are shown):
head(X) MouseType gender Caffeine RPM7
1 Control runner Male Red Bull 10.5461
2 High runner Female Gatorade 24.2718
3 High runner Female Gatorade 19.7215
4 Control runner Male Red Bull 13.5099
5 Control runner Male Red Bull 15.3210
6 High runner Male Gatorade 30.7953To get the experimental design, the table() function nicely summarize the numbers within each group:
table(X$gender:X$MouseType, X$Caffeine)
Water Gatorade Red Bull
Male:Control runner 10 14 12
Male:High runner 29 31 29
Female:Control runner 11 13 9
Female:High runner 28 32 31Before doing anything else, the data is plotted to visualize the response as function of the design (figure 8.1).
ggplot(X, aes(Caffeine, RPM7, color = Caffeine)) +
geom_boxplot() +
geom_jitter(alpha = .75) + # Alpha controls opacity
facet_wrap(~MouseType) +
theme_bw() +
theme(legend.position = "none") # Remove the legend
As an example, the mean and standard deviation is calculated for females of the unselected breed given water. The example is given using both base R and Tidyverse.
First, the information is extracted from the dataset:
extracted_x <- X[
X$MouseType == "Control runner" & X$gender == "Female" & X$Caffeine == "Water", # Rows
"RPM7"] # ColumnAnd then the extracted information can be fed into the mean() and sd() functions.
mean(extracted_x)[1] 8.649382sd(extracted_x)[1] 2.07804In Tidyverse the filter-function is used to filter by the MouseType, gender, and Caffeine levels of interest. Afterwards, the mean and standard deviation are computed using summarise.
X |>
filter(
MouseType == "Control runner",
gender == "Female",
Caffeine == "Water"
) |>
summarise(
mean = mean(RPM7),
sd = sd(RPM7)
) mean sd
1 8.649382 2.07804Both methods output a mean and standard deviation of
\[ \bar{x} = \sum_{i=1}^{n} \frac{x_i}{n} \approx 8.6 \]
\[ \hat{\sigma} = \sum_{i=1}^{n} \frac{(x_i - \bar{x})}{n - 1} \approx 2.1 \]
for the particular group of interest.
A more efficient way of doing this would be to calculate the mean and standard deviation of each combination and display them in a table. The example below is given using both base R and Tidyverse.
In base R aggregate() can be used to calculate a function over a group.
# Compute number of observations, mean and sd across groups
n <- aggregate(X$RPM7, by = list(X$MouseType, X$gender, X$Caffeine), length)
m <- aggregate(X$RPM7, by = list(X$MouseType, X$gender, X$Caffeine), mean)
s <- aggregate(X$RPM7, by = list(X$MouseType, X$gender, X$Caffeine), sd)
# Combine and set column names
stat_table <- cbind(n, m$x, s$x)
colnames(stat_table) <- c("MouseType","Gender","Caffeine","N","mean","sd")
# Print as nice markdown table
kable(stat_table, digits = 2)| MouseType | Gender | Caffeine | N | mean | sd |
|---|---|---|---|---|---|
| Control runner | Male | Water | 10 | 8.55 | 1.59 |
| High runner | Male | Water | 29 | 26.03 | 6.35 |
| Control runner | Female | Water | 11 | 8.65 | 2.08 |
| High runner | Female | Water | 28 | 24.59 | 7.25 |
| Control runner | Male | Gatorade | 14 | 8.55 | 2.09 |
| High runner | Male | Gatorade | 31 | 25.04 | 5.52 |
| Control runner | Female | Gatorade | 13 | 9.31 | 2.42 |
| High runner | Female | Gatorade | 32 | 22.65 | 5.86 |
| Control runner | Male | Red Bull | 12 | 10.73 | 2.95 |
| High runner | Male | Red Bull | 29 | 24.81 | 3.65 |
| Control runner | Female | Red Bull | 9 | 10.18 | 2.14 |
| High runner | Female | Red Bull | 31 | 24.50 | 6.50 |
In Tidyverse, the group_by()-function is used to group the observations by everything but the column RPM7 (so that will be type, gender and caffeine). Then, the number of observations, mean and standard devations are calculated across groups, giving the information for each group in one go.
X |>
group_by(across(-RPM7)) |>
summarise(
N = n(),
mean = mean(RPM7),
sd = sd(RPM7)
) |>
kable(digits = 2) # Output pretty table in markdown`summarise()` has grouped output by 'MouseType', 'gender'. You can override
using the `.groups` argument.| MouseType | gender | Caffeine | N | mean | sd |
|---|---|---|---|---|---|
| Control runner | Male | Water | 10 | 8.55 | 1.59 |
| Control runner | Male | Gatorade | 14 | 8.55 | 2.09 |
| Control runner | Male | Red Bull | 12 | 10.73 | 2.95 |
| Control runner | Female | Water | 11 | 8.65 | 2.08 |
| Control runner | Female | Gatorade | 13 | 9.31 | 2.42 |
| Control runner | Female | Red Bull | 9 | 10.18 | 2.14 |
| High runner | Male | Water | 29 | 26.03 | 6.35 |
| High runner | Male | Gatorade | 31 | 25.04 | 5.52 |
| High runner | Male | Red Bull | 29 | 24.81 | 3.65 |
| High runner | Female | Water | 28 | 24.59 | 7.25 |
| High runner | Female | Gatorade | 32 | 22.65 | 5.86 |
| High runner | Female | Red Bull | 31 | 24.50 | 6.50 |
Example Example 8.1 estimates the activity distribution of mice based on experimental data. Under the assumption that these metrics, the mean and the standard deviation, perfectly describe the distribution of the data, and further that this distribution is the normal distribution, we wish to calculate the probability of a given female mice, from normal breed, assigned water to perform higher than 10 RPM7.
That is:
\[\begin{align*} P(X \geq 10) &= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{10}^{\infty} e^{-\frac{(z - \mu)^2}{2\sigma}} dz \\ &= 1 - P(X < 10) \\ %X &= X \\ %X &= X &= 1 - \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^{10} e^{-\frac{(z - \mu)^2}{2\sigma}} dz \\ %N(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-\frac{1}{2}z^2} dz \end{align*}\]
where \(X\) is assummed normally distributed with mean \(\mu\) and standard deviation \(\sigma\).
That is: \(X \sim \mathcal{N} (\mu,\sigma^2)\). \
A small illustration below (figure 8.2) shows the area of interest.
# Create normal distributed data using estimated mean and sd
df <- data.frame(x = seq(2, 16, 0.01)) |>
mutate(y = dnorm(x, mean(extracted_x), sd(extracted_x)))
# Plot normal curve
ggplot(df, aes(x, y)) +
geom_line(linewidth = 1) +
geom_area(data = subset(df, x >= 10), fill = "red") +
labs(
x = "RPM7",
y = "Probability"
) +
theme_bw()
The size of the red area is simply calculated via the pnorm() function.
1 - pnorm(10, mean(extracted_x), sd(extracted_x))[1] 0.2578628\[ P(X \geq 10) \approx 0.26 \]
The probability of a single female mice, from normal breed, administered with water to perform higher than \(10 RPM7\) is hence \(0.26\).
The confidence interval for the center of a normal distribution (\(\mu\)) is calculated as follows:
\[ CI_{\mu,1-\alpha}: \hat{\mu} \pm t_{1-\alpha/2,df} \cdot \hat{\sigma} / \sqrt(n) \] \[ \bar{X} \pm t_{1-\alpha/2,df} \cdot s / \sqrt(n) \]
where \(t_{1-\alpha/2,df}\) is a fractile, a number, which determines the coverage. Here, \(\alpha\) is the left ot part. I.e. if a \(90 \%\) confidence interval is wanted, the left out part is \(\alpha = 0.10\). \(n\) is the number of samples on which the mean (\(\bar{X}\)) is estimated, and \(df\) is the degrees of freedom, which refers to how well the standard deviation is estimated.
For instance, if one needs a \(95 \%\) confidence interval (\(\alpha = 0.05\)) based on a sample of \(20\) observations \(t_{1-\alpha/2,df} = t_{0.975,20-1} = 2.093\).
It is further noticed that the interval is symmetric around the mean, and that the more samples (\(n\)) the lower the spread as the standard deviation is divided by \(\sqrt{n}\).
In example Example 8.1 concerning the relation between caffeine and voluntary activity of mice, the running activity of mice in a wheel is recorded. In this example we wish to estimate the mean activity of the two groups with either water- or Red bull as caffeine source in control runner female mice, and give a confidence interval for these means.
First, the data is subset to only include the two groups of interest.
x_subset <- X |>
filter(
MouseType == "Control runner",
gender == "Female",
Caffeine %in% c("Red Bull", "Water")
)After subsetting the data to only include these two groups a plot of the raw data looks like the following, and consists of \(n = n_1 + n_2 = 11+ 9 = 20\) samples (figure 8.3).
ggplot(x_subset, aes(Caffeine, RPM7, color = Caffeine)) +
geom_boxplot() +
geom_jitter() +
theme_bw() +
theme(legend.position = "none") # Remove legend
The mean and standard deviation for the two groups can be calculated using summarise() via Tidyverse or aggregate() in base R (see Example 8.1 for an example on aggregate).
x_subset |>
group_by(Caffeine) |>
summarise(
N = n(),
mean = mean(RPM7),
sd = sd(RPM7)
) |>
kable(digits = 2) # For pretty tables in markdown| Caffeine | N | mean | sd |
|---|---|---|---|
| Water | 11 | 8.65 | 2.08 |
| Red Bull | 9 | 10.18 | 2.14 |
The confidence interval is calculated as follows:
\[\begin{align} CI_{\mu}: & \hat{\mu} \pm t_{0.975,n-1} \cdot \hat{\sigma} / \sqrt(n) \\ & \bar{X} \pm t_{0.975,n-1} \cdot s / \sqrt(n) \end{align}\]
With a total of \(df_1 = n_1 - 1 = 10\) degress of freedom the t-fractileat level \(95 \%\) is \(t_{0.975,10} = 2.23\) why the confidence interval for the water treated group is:
\[ CI_{\mu_{water}}: [8.6 - 2.23 \cdot 2.1 / \sqrt{11}; 8.6 + 2.23 \cdot 2.1/ \sqrt{11}] = [7.3; 10.0] \]
n <- 11 # Number of obs in water group
t_frac <- qt(1 - 0.05/2, n - 1) # T-value
s <- sd(x_subset[x_subset$Caffeine == "Water", "RPM7"])
m <- mean(x_subset[x_subset$Caffeine == "Water", "RPM7"])
ci_lower <- m - t_frac * s / sqrt(n)
ci_upper <- m + t_frac * s / sqrt(n)
cat("95% confidence interval for water group:\n", c(ci_lower, ci_upper))95% confidence interval for water group:
7.253336 10.04543And similar for the Red bull treated group.
\[\begin{equation} CI_{\mu_{Red Bull}}: [8.5; 11.8] \end{equation}\]
In R this can be assesed via the t.test() function.
t.test(x_subset[x_subset$Caffeine == "Water", "RPM7"])
One Sample t-test
data: x_subset[x_subset$Caffeine == "Water", "RPM7"]
t = 13.805, df = 10, p-value = 7.744e-08
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
7.253336 10.045428
sample estimates:
mean of x
8.649382 Height and weight are important measures of growth during childhood. However, the largest factor impacting these measures are naturally the age of the child. In order to be able to compare children with slightly different ages such data are transformed using so-called growth curves. These are provided by the WHO and are constructed to represent the world wide distribution at specific ages for boys and girl. In this exercise you are supposed to take height and weight measures from African children living under different circumstances, transform them using WHO numbers for the world wide distribution, and further look for explanatory variables causing differences in growth.
You can use the command X <- X[complete.cases(X),] to remove the empty rows in the dataframe X, or google how to compute descriptive statistics in the case of missing values.
This can be done by creating a new vector inserting the length-vector from the previous dataframe to the equation for Z-score.