10.4: May the (Strong) Force be with You

Ruritania is a patron of high-energy physics — and of Star Wars enthusiasts. King Rudolph donated several million crowns to Switzerland to aid in researching the strong force. That money was used at CERN (the Conseil Européen pour la Recherche Nucléaire) for several experiments. Each experiment consisted of a beam of protons crashing into a target. That beam had a constant energy level. What changed were the target sizes and the energy level of the proton after the collision. Many experiments were run at each energy level, and the standard deviation of the energies was measured.

In a theory proposed by Ruritanian scientists that is not entirely clear to His Majesty (or to your author), there should be a linear relationship between the cross sectional area and the inverse of the energy. The data are given in the table below.

Table \(\PageIndex{1}\): The data table for this example.

Cross Section[b]	Energy [MeV]	StDev [MeV]
1	848.9	7.8
2	476.9	9.2
3	350.9	9.4
4	289.2	10.2
5	251.7	7.4
6	225.8	9.3
7	209.79	7.2
8	193.9	5.3

Table \(\PageIndex{1}\): Data for the example regarding the strong nuclear force. Units are given in brackets.

The first column is the value of the independent variable. The second column is the mean of the energy level of the photon after the collision. The third column is the standard deviation in those energy levels. Note that the variability at each cross-section differs. This is based on both the number of experiments and the inherent variability at that area.

Ordinary Least Squares

Let us ignore the different uncertainties in each energy level (the standard deviations). That is, let us just fit this as an OLS model.

Here is the code:

barns  = c(1,2,3,4,5,6,7,8)
energy = c(848.9,476.9,350.9,289.2,251.7,225.8,209.7,193.9)
Ibarns = 1/barns

modOLS = lm(energy~Ibarns)

summary(modOLS)
confint(modOLS)

The output suggests that the relationship between the cross sectional area of the target and the inverse of the resulting energy of the photon is statistically significant.

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 101.9148     0.5464   186.5 1.60e-12 ***
Ibarns      747.5306     1.2505   597.8 1.48e-15 ***
---

Residual standard error: 0.9719 on 6 degrees of freedom
Multiple R-squared:      1,    Adjusted R-squared:      1
F-statistic: 3.574e+05 on 1 and 6 DF,  p-value: 1.479e-15

A 95% confidence interval for the relationship is from 744.5 to 750.6 (with units of MeV\(\cdot\)barn).

Weighted Least Squares

Note, however, that the uncertainty in the measurements varies. We are more uncertain with some of our estimated energy than with others. If we do not take this uncertainty into consideration, we may be biasing our results. To include this information, we can use weighted least squares regression. The code to fit with weighted least regression is as follows:

barns  = c(1,2,3,4,5,6,7,8)
energy = c(848.9,476.9,350.9,289.2,251.7,225.8,209.7,193.9)
stdev  = c(7.8,9.2,9.4,10.2,7.4,9.3,7.2,5.3)
v      = stdev^2

Ibarns = 1/barns

modWLS = lm(energy~Ibarns, weight=1/v)
summary(modWLS)
confint(modWLS)

Note that we are 95% confident the effect of the target's cross section on the resulting energy is from 744.5 to 751.2 MeV\(\cdot\)barn.

Important R Feature

In R, as in many statistical programs, the weights you provide in the function call are inversely proportional to the variances. This is why we used weight=1/v in the function call.

With that being said, always check the documentation to make doubly sure. Frequently, this information is difficult to find and an error will not be thrown to let you know.

How do we know this?

Or, an even better question:

How do we check that the weights really are proportional to the inverse of the variance?

The short answer is to check your work "by hand." As we once did OLS by hand using the matrix functions in R, we can also do WLS by hand. As always, make sure you know what each line does and why they are put together as they are:

## The data
barns  = c(1,2,3,4,5,6,7,8)
energy = c(848.9,476.9,350.9,289.2,251.7,225.8,209.7,193.9)
stdev  = c(7.8,9.2,9.4,10.2,7.4,9.3,7.2,5.3)

## A few minor calculations
v = stdev^2
Ibarns = 1/barns
n = length(Ibarns)

## The needed matrices
Y = matrix(energy, ncol=1)
X = matrix( c(rep(1,n), Ibarns), ncol=2 )
colnames(X) = c("b0","b1")
D = diag(v)

## The estimate vector
solve(t(X) %*% solve(D) %*% X) %*% t(X) %*% solve(D) %*% Y

When I run this, I get the following output

       [,1]
b0 101.6074
b1 747.8364

These are the matrix calculations according to the WLS estimation formula. Thus, this is the correct answer. Double-check that this known correct answer matches the answer given in modWLS. If it does not, then you will need to change the expression in the weights part of the function call. In R, it will match. In other pieces of software, it may not.

Be aware!

Caution

The difference between the effects estimated from using ordinary least squares and using weighted least squares is rather minor in this example. It need not be, as the next example shows.