There are four exercises below. Three exercises are worth 10 points each, and are labelled (Either R, Python or Julia). For these exercises, you can implement a solution in the language of your choice. I’ve included code chunk templates for all three languages, but you only need to fill in the chunks for one language.
One exercise is labelled (R, Python and Julia). This exercise is worth 20 points, and you must provide a solution using all three languages. This exercise also includes code chunk templates for all three languages.
Some of the exercises have questions relating to the results. Write your answers in the narrative part of the text, formatted as italicized text, bold text or
quoted text.
Submit both Rmd and typeset results using file name
formats as in previous exercises.
Do not use libraries in R that are not part of the base R
installation; this includes popular libraries for data tables like
data.table or tidyverse libraries like dplyr
and tidyr. You may use any Python or Julia libraries
introduced in lecture, except for graphs. You are free to use any Python
or Julia graphics library, as long as I can install it on my machine if
necessary.
The OzDASL web site includes a data set titled Pain Thresholds of Blonds and Brunettes (http://www.statsci.org/data/oz/blonds.html). The data are as follows:
| HairColour | Pain |
|---|---|
| LightBlond | 62 |
| LightBlond | 60 |
| LightBlond | 71 |
| LightBlond | 55 |
| LightBlond | 48 |
| DarkBlond | 63 |
| DarkBlond | 57 |
| DarkBlond | 52 |
| DarkBlond | 41 |
| DarkBlond | 43 |
| LightBrunette | 42 |
| LightBrunette | 50 |
| LightBrunette | 41 |
| LightBrunette | 37 |
| DarkBrunette | 32 |
| DarkBrunette | 39 |
| DarkBrunette | 51 |
| DarkBrunette | 30 |
| DarkBrunette | 35 |
Copy the data from the table above to create a
data.frame in R, a pandas
dataframe in Python, or a DataFrame in Julia.
Write the code the appropriate chunk below. Name the data table
PainThreshold. Note that the values in column one will need
to be coded as text (with single or double quotes as necessary). Don’t
download the linked text file; this is an exercise in programmatically
creating a data table. Take care to format the data table code so that
it fits in the typeset document.
After creating the table, modify HairColour so that it
is an ordinal variable, with the order LightBlond,
DarkBlond, LightBrunette and
DarkBrunette.
# Creating table
HairColour <- c("LightBlond", "LightBlond", "LightBlond", "LightBlond", "LightBlond",
"DarkBlond", "DarkBlond", "DarkBlond", "DarkBlond", "DarkBlond",
"LightBrunette", "LightBrunette", "LightBrunette", "LightBrunette",
"DarkBrunette", "DarkBrunette", "DarkBrunette", "DarkBrunette", "DarkBrunette")
Pain <- c(62, 60, 71, 55, 48, 63, 57, 52, 41, 43, 42, 50, 41, 37, 32, 39, 51, 30, 35)
PainThreshold <- data.frame(HairColour, Pain)
# Modifying HairColour to be an ordinal variable
PainThreshold$HairColour <- factor(PainThreshold$HairColour,
levels = c("LightBlond", "DarkBlond", "LightBrunette", "DarkBrunette"),
ordered = TRUE)
print(PainThreshold)
## HairColour Pain
## 1 LightBlond 62
## 2 LightBlond 60
## 3 LightBlond 71
## 4 LightBlond 55
## 5 LightBlond 48
## 6 DarkBlond 63
## 7 DarkBlond 57
## 8 DarkBlond 52
## 9 DarkBlond 41
## 10 DarkBlond 43
## 11 LightBrunette 42
## 12 LightBrunette 50
## 13 LightBrunette 41
## 14 LightBrunette 37
## 15 DarkBrunette 32
## 16 DarkBrunette 39
## 17 DarkBrunette 51
## 18 DarkBrunette 30
## 19 DarkBrunette 35
The linked source for these data includes a box-whisker plot, with
HairColour as the independent variable. Reproduce the
box-whisker plot below. Did you enter the data correctly? Compare with
the source plot.
The visuals look very similar, these differences are subjective to what the user wants (colors, size, etc). Overall the data shares the same analysis of the source plot.
“Pain threshold decreases as hair darkness increases (blonds are tougher!). Light blonds are significantly different from both brunette categories. Other contrasts are not significant.”
# Create box-whisker plot using base R
boxplot(Pain ~ HairColour, data = PainThreshold,
main = "Pain Thresholds by Hair Colour",
xlab = "Hair Colour", ylab = "Pain Threshold",
col = "white", border = "black")
You will be asked to use your ConfidenceInterval
function for this exercise. Include the function definition here.
# function definitions
def ConfidenceInterval(data, confidence=0.95):
n = len(data)
mean = np.mean(data)
stderr = np.std(data, ddof=1) / np.sqrt(n)
error = stats.t.ppf((1 + confidence) / 2., n - 1) * stderr
return (mean - error, mean + error)
# function definitions
Go to http://www.itl.nist.gov/div898/strd/anova/SiRstv.html
and use the data listed under Data File in Table Format (https://www.itl.nist.gov/div898/strd/anova/SiRstvt.dat)
Edit this into a file (tab delimited, .csv, etc,) that
can be read into R, Python or Julia, or find an appropriate function
that can read the file as-is. You will need to upload the edited file to
D2L along with your Rmd files. Provide a brief comment on changes you
make, or assumptions about the file needed for you file to be read into
R, Python or Julia. Read the data into a data table.
# Define the URL and read file into R
url <- "https://www.itl.nist.gov/div898/strd/anova/SiRstvt.dat"
data <- read.table(url, header = FALSE, sep = "", fill = TRUE)
# Path to save CSV
save_path <- "C:/Users/Allen/OneDrive - Dakota State University/Summer 24/Statistical Programming 600/Week 5/HW/SiRstvt.csv"
# Save the data to the specified path
write.csv(data, file = save_path, row.names = FALSE)
# Print data
print(data)
## V1 V2 V3 V4
## 1 NIST/ITL StRD
## 2 Dataset Name: SiRstv (SiRstvt.dat)
## 3 File Format: ASCII
## 4 Certified Values (lines 41
## 5 Data (lines 61 to
## 6 Procedure: Analysis of Variance
## 7 Reference: Ehrstein, James and
## 8 Carroll.
## 9 Unpublished NIST dataset.
## 10 Data: 1 Factor
## 11 5 Treatments
## 12 5 Replicates/Cell
## 13 25 Observations
## 14 3 Constant Leading Digits
## 15 Lower Level of Difficulty
## 16 Observed Data
## 17 Model: 6 Parameters (mu,tau_1,
## 18 tau_5)
## 19 y_{ij} = mu +
## 20 epsilon_{ij}
## 21 Certified Values:
## 22 Source of Sums of
## 23 Variation df Squares Squares
## 24 Between Instrument 4 5.11462616000000E-02
## 25 Within Instrument 20 2.16636560000000E-01
## 26 Certified R-Squared 1.90999039051129E-01
## 27 Certified Residual
## 28 Standard Deviation 1.04076068334656E-01
## 29 Data:
## 30 Instrument
## 31 1 2 3 4
## 32 196.3052 196.3042 196.1303 196.2795
## 33 196.1240 196.3825 196.2005 196.1748
## 34 196.1890 196.1669 196.2889 196.1494
## 35 196.2569 196.3257 196.0343 196.1485
## 36 196.3403 196.0422 196.1811 195.9885
## V5 V6
## 1
## 2
## 3
## 4 to 47)
## 5 65)
## 6
## 7 Croarkin, M.
## 8
## 9
## 10
## 11
## 12
## 13
## 14
## 15
## 16
## 17 ... ,
## 18
## 19 tau_i +
## 20
## 21
## 22 Mean
## 23 F Statistic
## 24 1.27865654000000E-02 1.18046237440255E+00
## 25 1.08318280000000E-02
## 26
## 27
## 28
## 29
## 30
## 31 5
## 32 196.2119
## 33 196.1051
## 34 196.1850
## 35 196.0052
## 36 196.2090
I assume that the data file at the provided URL
("https://www.itl.nist.gov/div898/strd/anova/SiRstvt.dat")
is tab-delimited. The first row appears to have header, but issues
persist to remove the header;(header = TRUE) would initiate
the write functionality. With that said I managed to initiate the CSV by
setting the header to “FALSE” and can later remove the headers
manually.
There are 5 columns in these data. Calculate mean, standard deviation and sample size for each column in this data, using column summary functions. Print the results below.
import numpy as np
import scipy.stats as stats
import pandas as pd
from pprint import pprint
# Specify the path to your CSV file
file_path = r"C:\Users\Allen\OneDrive - Dakota State University\Summer 24\Statistical Programming 600\Week 5\HW\SiRstvt.csv"
# Read the CSV file into a DataFrame
data = pd.read_csv(file_path)
# Extracting relevant rows and columns
cleaned_data = data.iloc[31:35, :5]
cleaned_data.columns = ['V1', 'V2', 'V3', 'V4', 'V5']
# Convert columns to numeric, forcing errors to NaN (which will be handled later)
cleaned_data = cleaned_data.apply(pd.to_numeric, errors='coerce')
# Display cleaned data
print("Table:")
## Table:
print(cleaned_data)
## V1 V2 V3 V4 V5
## 31 196.3052 196.3042 196.1303 196.2795 196.2119
## 32 196.1240 196.3825 196.2005 196.1748 196.1051
## 33 196.1890 196.1669 196.2889 196.1494 196.1850
## 34 196.2569 196.3257 196.0343 196.1485 196.0052
# Calculate statistics for each column
summary_stats = cleaned_data.apply(lambda col: {
'mean': np.mean(col),
'sd': np.std(col, ddof=1),
'n': len(col[~np.isnan(col)])
})
# Method to print the full view of the stats!
print("\nSummary Statistics:")
##
## Summary Statistics:
pprint(summary_stats.to_dict())
## {'V1': {'mean': 196.218775, 'n': 4, 'sd': 0.07914469344183006},
## 'V2': {'mean': 196.29482499999997, 'n': 4, 'sd': 0.09145648783255737},
## 'V3': {'mean': 196.1635, 'n': 4, 'sd': 0.10784099406070445},
## 'V4': {'mean': 196.18805, 'n': 4, 'sd': 0.06217365465640073},
## 'V5': {'mean': 196.1268, 'n': 4, 'sd': 0.09289187262619278}}
Reuse your ConfidenceInterval function to compute
confidence intervals for the means in this data set. Note, you can do
this with one function call if you use vectors, list comprehensions or
broadcasting.
# Compute confidence intervals for each column
confidence_intervals = cleaned_data.apply(ConfidenceInterval)
# Pretty-print the confidence intervals
print("\nConfidence Intervals:")
##
## Confidence Intervals:
pprint(confidence_intervals.to_dict())
## {'V1': {0: 196.0928381313866, 1: 196.3447118686134},
## 'V2': {0: 196.14929731910146, 1: 196.44035268089849},
## 'V3': {0: 195.99190091344667, 1: 196.33509908655333},
## 'V4': {0: 196.08911784122637, 1: 196.28698215877364},
## 'V5': {0: 195.97898830158493, 1: 196.27461169841507}}
We will use data from https://acsess.onlinelibrary.wiley.com/doi/abs/10.2134/jeq2007.0099, Table 1. The original paper is also available on D2L.
Download the file Khan.csv from D2L and read the file
into a data frame. Print a summary of the table.
import pandas as pd
# Path to your CSV file
file_path = r"C:\Users\Allen\OneDrive - Dakota State University\Summer 24\Statistical Programming 600\Week 5\HW\Khan.csv"
# Read the CSV file into a df
df = pd.read_csv(file_path)
# Display the first few rows of the df
print(df.head())
## Rotation Fertilizer Depth Mean55 Mean05 SD05 Diff
## 0 C-C none 0-15 1.376 1.168 0.007 -0.208
## 1 C-C none 15-30 1.342 1.068 0.008 -0.274
## 2 C-C none 30-46 1.020 0.996 0.006 -0.024
## 3 C-C NPK 0-15 1.376 1.268 0.008 -0.108
## 4 C-C NPK 15-30 1.342 1.210 0.010 -0.132
To show that the data was read correctly, create three plots. Plot
Mean05 and Mean55 are the amount of soil
organic carbon measured in crop land experimental units in 2005 and 1955
respectively. Rotation is the crop rotation plan (i.e. corn
followed by soybeans followed by corn) for the respective plots, and
Fertilizer is the type of fertilizer applied to the plots
over the period from 1955-2005.
These three plots should reproduce the three types of plots shown in
the Week 5 Data Tables Plots video, Categorical vs
Categorical, Continuous vs Continuous and
Continuous vs Categorical. Add these as titles to your
plots, as appropriate. If you choose Julia for this exercise, you will
not be required to plot Categorical vs Categorical.
Do you notice anything unusual about the data?
import pandas as pd
import matplotlib.pyplot as plt
# 1. Rotation vs Fertilizer (Categorical vs Categorical)
plt.figure(figsize=(12, 6))
rotation_fertilizer_ct = pd.crosstab(df['Rotation'], df['Fertilizer'])
rotation_fertilizer_ct.plot(kind='bar', stacked=True, colormap='viridis')
plt.title('Categorical vs Categorical: Rotation vs Fertilizer')
plt.xlabel('Rotation')
plt.ylabel('Count')
plt.xticks(rotation=45)
## (array([0, 1, 2]), [Text(0, 0, 'C-C'), Text(1, 0, 'C-O(S)'), Text(2, 0, 'C-O-H')])
plt.legend(title='Fertilizer')
plt.tight_layout()
plt.show()
# 2. Mean55 vs Fertilizer (Continuous vs Categorical)
plt.figure(figsize=(12, 6))
ax = df.boxplot(column='Mean55', by='Fertilizer', grid=False)
ax.set_title('Continuous vs Categorical: Mean55 vs Fertilizer')
plt.suptitle('') # Suppress the default title to make it cleaner
ax.set_xlabel('Fertilizer')
ax.set_ylabel('Mean55')
plt.tight_layout()
plt.show()
# 3. Mean55 vs Mean05 (Continuous vs Continuous)
plt.figure(figsize=(12, 6))
plt.scatter(df['Mean05'], df['Mean55'], alpha=0.5)
plt.title('Continuous vs Continuous: Mean55 vs Mean05')
plt.xlabel('Mean05')
plt.ylabel('Mean55')
plt.grid(True)
plt.tight_layout()
plt.show()
1. Rotation vs Fertilizer (Categorical vs Categorical): Stacked bar
plot; Each rotation type (C-C, C-O(S), C-O(H)) has a consistent
distribution of the three types of fertilizers (HNPK, NPK, none). There
is no significant variation in the count of fertilizers across different
rotations.
2. Mean55 vs Fertilizer (Continuous vs Categorical): Box plot; The distribution of Mean55 values varies with the type of fertilizer used. HNPK fertilizer shows the highest median and interquartile range. NPK and none fertilizers show lower medians and more similar distributions. There are no extreme outliers, but there is noticeable variation among the groups.
3. Mean55 vs Mean05 (Continuous vs Continuous): Scatter plot; There appears to be a positive correlation between Mean05 and Mean55. As the Mean05 value increases, the Mean55 value tends to increase as well. However, there are some points where the Mean05 value is low (around 1.0) but Mean55 values are higher, indicating potential outliers or errors in the data.
Example data for the chemical composition of glass fragments are
available https://archive.ics.uci.edu/ml/datasets/Glass+Identification.
I’ve uploaded the data to D2L in the file glass.data. This
is in CSV format with no header. Read the data into R, Python or
Julia.
Add the following column headers to the data table (see the source link for a description of the data columns).
IdRINaMgAlSiKCaBaFeTypeYou should not print the data tables in the typeset document. Instead, steps in Part b will tell us if the columns have been properly renamed. However, you should print the number of rows in the data. There should be 214 rows and 11 columns.
# Define the path to the file
file_path <- "C:/Users/Allen/OneDrive - Dakota State University/Summer 24/Statistical Programming 600/Week 5/HW/glass.data"
# Define column names
columns <- c('Id', 'RI', 'Na', 'Mg', 'Al', 'Si', 'K', 'Ca', 'Ba', 'Fe', 'Type')
# Read the data into a data frame
glass_df <- read.csv(file_path, header=FALSE, col.names=columns)
# Print the number of rows and columns
cat("Number of rows:", nrow(glass_df), "Number of columns:", ncol(glass_df), "\n")
## Number of rows: 214 Number of columns: 11
import pandas as pd
# Define the path to the file
file_path = r"C:\Users\Allen\OneDrive - Dakota State University\Summer 24\Statistical Programming 600\Week 5\HW\glass.data"
# Define column names
columns = ['Id', 'RI', 'Na', 'Mg', 'Al', 'Si', 'K', 'Ca', 'Ba', 'Fe', 'Type']
# Read the data into a DataFrame
glass_df = pd.read_csv(file_path, header=None, names=columns)
# Print the number of rows and columns
print(f'Number of rows: {glass_df.shape[0]}, Number of columns: {glass_df.shape[1]}')
## Number of rows: 214, Number of columns: 11
using DataFrames
using CSV
# Define the path to the file
file_path = "C:/Users/Allen/OneDrive - Dakota State University/Summer 24/Statistical Programming 600/Week 5/HW/glass.data"
## "C:/Users/Allen/OneDrive - Dakota State University/Summer 24/Statistical Programming 600/Week 5/HW/glass.data"
# Define column names
columns = [:Id, :RI, :Na, :Mg, :Al, :Si, :K, :Ca, :Ba, :Fe, :Type]
## 11-element Vector{Symbol}:
## :Id
## :RI
## :Na
## :Mg
## :Al
## :Si
## :K
## :Ca
## :Ba
## :Fe
## :Type
# Read the data into a DataFrame
glass_df = CSV.read(file_path, DataFrame; header=false)
## 214×11 DataFrame
## Row │ Column1 Column2 Column3 Column4 Column5 Column6 Column7 Column8 ⋯
## │ Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 ⋯
## ─────┼──────────────────────────────────────────────────────────────────────────
## 1 │ 1 1.52101 13.64 4.49 1.1 71.78 0.06 8.75 ⋯
## 2 │ 2 1.51761 13.89 3.6 1.36 72.73 0.48 7.83
## 3 │ 3 1.51618 13.53 3.55 1.54 72.99 0.39 7.78
## 4 │ 4 1.51766 13.21 3.69 1.29 72.61 0.57 8.22
## 5 │ 5 1.51742 13.27 3.62 1.24 73.08 0.55 8.07 ⋯
## 6 │ 6 1.51596 12.79 3.61 1.62 72.97 0.64 8.07
## 7 │ 7 1.51743 13.3 3.6 1.14 73.09 0.58 8.17
## 8 │ 8 1.51756 13.15 3.61 1.05 73.24 0.57 8.24
## ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
## 208 │ 208 1.51831 14.39 0.0 1.82 72.86 1.41 6.47 ⋯
## 209 │ 209 1.5164 14.37 0.0 2.74 72.85 0.0 9.45
## 210 │ 210 1.51623 14.14 0.0 2.88 72.61 0.08 9.18
## 211 │ 211 1.51685 14.92 0.0 1.99 73.06 0.0 8.4
## 212 │ 212 1.52065 14.36 0.0 2.02 73.42 0.0 8.44 ⋯
## 213 │ 213 1.51651 14.38 0.0 1.94 73.61 0.0 8.48
## 214 │ 214 1.51711 14.23 0.0 2.08 73.36 0.0 8.62
## 3 columns and 199 rows omitted
rename!(glass_df, columns)
## 214×11 DataFrame
## Row │ Id RI Na Mg Al Si K Ca B ⋯
## │ Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 F ⋯
## ─────┼──────────────────────────────────────────────────────────────────────────
## 1 │ 1 1.52101 13.64 4.49 1.1 71.78 0.06 8.75 ⋯
## 2 │ 2 1.51761 13.89 3.6 1.36 72.73 0.48 7.83
## 3 │ 3 1.51618 13.53 3.55 1.54 72.99 0.39 7.78
## 4 │ 4 1.51766 13.21 3.69 1.29 72.61 0.57 8.22
## 5 │ 5 1.51742 13.27 3.62 1.24 73.08 0.55 8.07 ⋯
## 6 │ 6 1.51596 12.79 3.61 1.62 72.97 0.64 8.07
## 7 │ 7 1.51743 13.3 3.6 1.14 73.09 0.58 8.17
## 8 │ 8 1.51756 13.15 3.61 1.05 73.24 0.57 8.24
## ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
## 208 │ 208 1.51831 14.39 0.0 1.82 72.86 1.41 6.47 ⋯
## 209 │ 209 1.5164 14.37 0.0 2.74 72.85 0.0 9.45
## 210 │ 210 1.51623 14.14 0.0 2.88 72.61 0.08 9.18
## 211 │ 211 1.51685 14.92 0.0 1.99 73.06 0.0 8.4
## 212 │ 212 1.52065 14.36 0.0 2.02 73.42 0.0 8.44 ⋯
## 213 │ 213 1.51651 14.38 0.0 1.94 73.61 0.0 8.48
## 214 │ 214 1.51711 14.23 0.0 2.08 73.36 0.0 8.62
## 3 columns and 199 rows omitted
# Print the number of rows and columns
println("Number of rows: $(size(glass_df, 1)), Number of columns: $(size(glass_df, 2))")
## Number of rows: 214, Number of columns: 11
Create two new data columns for \(\log(Ca/K)\) and \(\log(Ca/Si)\). The first is the log of
(column Ca divided by column K), the second is
the log of (column Ca divided by column Si).
You don’t need to use the column names log(Ca/K) or
log(Ca/Si), you can use names without special characters
(i.e. ‘(’ or ‘/’)). You should not print the data tables in the typeset
document. We’ll check the calculations in Part c.
# Calculate log(Ca/K) and log(Ca/Si)
glass_df$log_Ca_K <- log(glass_df$Ca / glass_df$K)
glass_df$log_Ca_Si <- log(glass_df$Ca / glass_df$Si)
import numpy as np
# Calculate log(Ca/K) and log(Ca/Si)
glass_df['log_Ca_K'] = np.log(glass_df['Ca'] / glass_df['K'])
glass_df['log_Ca_Si'] = np.log(glass_df['Ca'] / glass_df['Si'])
using DataFrames
using Statistics
# Calculate log(Ca/K) and log(Ca/Si)
glass_df.log_Ca_K = log.(glass_df.Ca ./ glass_df.K)
## 214-element Vector{Float64}:
## 4.98246441712956
## 2.791931685082912
## 2.993164878048745
## 2.6686891272216298
## 2.685990483037478
## 2.534440584910277
## 2.6451960843135836
## 2.671119262074922
## 2.6960740100554945
## 2.6903506240028094
## ⋮
## Inf
## Inf
## 1.5235864041227325
## Inf
## 4.7427558489406545
## Inf
## Inf
## Inf
## Inf
glass_df.log_Ca_Si = log.(glass_df.Ca ./ glass_df.Si)
## 214-element Vector{Float64}:
## -2.104552185349938
## -2.2287914441637375
## -2.2387661072732405
## -2.1785324443240666
## -2.2034012492279484
## -2.201894915495671
## -2.191222649620389
## -2.1847413758849297
## -2.16152109849753
## -2.162090739614273
## ⋮
## -2.1373981270239013
## -2.132096365512774
## -2.4213636830676086
## -2.0423877911355626
## -2.0680754487597564
## -2.1630493155178683
## -2.1632140694608224
## -2.1610704362668027
## -2.141293743272608
To show your work, plot log(Ca/Si) versus
log(Ca/K) (independent variable). The x- and y-axes of your
plot should be labelled log(Ca/K) and
log(Ca/Si).
# Plot log(Ca/Si) versus log(Ca/K)
plot(glass_df$log_Ca_K, glass_df$log_Ca_Si,
xlab = "log(Ca/K)", ylab = "log(Ca/Si)",
main = "Logarithm of Ca/K vs Ca/Si")
import matplotlib.pyplot as plt
# Plot log(Ca/Si) versus log(Ca/K)
plt.figure(figsize=(8, 6))
plt.scatter(glass_df['log_Ca_K'], glass_df['log_Ca_Si'])
plt.xlabel('log(Ca/K)')
plt.ylabel('log(Ca/Si)')
plt.title('Logarithm of Ca/K vs Ca/Si')
plt.grid(True)
plt.show()
using Plots
# Plot log(Ca/Si) versus log(Ca/K)
scatter(glass_df.log_Ca_K, glass_df.log_Ca_Si,
xlabel = "log(Ca/K)", ylabel = "log(Ca/Si)",
title = "Logarithm of Ca/K vs Ca/Si")
As a side note, Dr. Saunders has been using a plot similar to this on the homepage for STAT 601 and STAT 602.