A Tale of "*z* and *p*"

When I was a young man, I traveled extensively to exotic foreign lands.
One of the most curious of these was New Mexico --- where, like Wyoming,
you don't even need a passport for your visit, and the local currency is
the dollar. These conveniences count for a lot. Think about this if you
are considering going to Tahiti or Bora Bora.

I recently returned to New Mexico for a visit and, while sipping a
cerveza (this is a local brew somewhat like beer, only better) on a café
patio in Old Town Albuquerque (spell that, if you can) I spied my old
friend and teacher, the brilliant statistician John McBabble. I offered
to buy him a cerveza, which no one has *ever* known John to turn down,
unless someone else is offering Scotch.

Our conversation soon turned to cosmic questions in statistics. John had
been working on the probability that the word *Albuquerque* would ever
be typed by a million monkeys on laptop computers, and had concluded
that it was about as likely as a flying pig. I had trouble following his
calculations, probably because he wrote them on a cerveza-soaked
cocktail napkin with the stylus from his palm pilot. John's statistical
results are never refuted.

I told McBabble that when I visited New Mexico in my younger days, in
the mid 1950's, Albuquerque had been as dry as his glass. He said that
he could help me prove that things were indeed incredibly dry in 1955 if
I'd buy him a refill, which I did as fast as a monkey can type.

Whipping out his laptop, and never at a loss for important data to
analyze, John showed me a table of data that he downloaded from the
Western Regional Climate Center. It was the historic rainfall data for
Albuquerque. I looked the table over, and asked whether it might not be
hard to work with, as it was monthly data. Then he said some complete
gibberish, "Not to worry, I R an expert in R," and laughed "R R R".

"Whoa," I thought, "time to hold back on the cervezas."

I watched as McBabble fired up a software program named "R," which he
kept saying over and over was free. You can count on a Scot to notice
things like that. Working in R, John changed R's working directory to
the folder where the data table was stored as a text file. Looking at
the text table again, I noted that the first row, which contained the
names of the months, had one less item in it (12) than the rest of the
rows, each of which started with the year and then had the rainfall
totals for each month (13 items). He said, "Watch this!" and typed

pcp.tbl\<-read.table("ABQprecip.txt",header=T)

pcp.tbl

and I saw the data all neatly lined up by year and month. And then he
typed

pcp.df\<-as.data.frame(pcp.tbl)

pcp.df

and I saw exactly the same table.

"That didn't seem to do much," I said.

"Oh, but now the data are in a data frame, which is much handier," he
replied.

I was at a loss to see how a data frame was handier than a table,
because they looked the same to me.

Not stopping to explain, McBabble forged ahead, typing

pcp.df\$Annual\<-rowSums(pcp.df);pcp.df

Looking over his shoulder, I saw that the data frame now had one more
column, labeled *Annual*, and that this column contained the row totals,
which were of course annual precipitations.

I said, "Cool, let's plot the annual data."

"Not so fast," he responded. "If you look closely, you can see that a
couple of years of data are missing, so if we want to plot the data as a
time series we'll need to do an *xy* plot with years on the *x* axis."

"What's the problem, the years are there in the row labels," I noted.

"True, but they're just labels and not numbers, so we can't use them for
the *x* variable."

"Well that's pretty awkward," I said. "Must be cheap software."

"Free, not cheap," he replied. "And it's easy to fix." He typed on:

pcp.df\$Year\<-as.numeric(row.names(pcp.df)) ; pcp.df

and I saw that the data frame had yet another column, labeled *Year*,
and I assumed that this column now held numbers.

Rubbing his hands together, and then turning back to the keyboard, John
typed

plot(pcp.df\$Year,pcp.df\$Annual,type="o")

and a graph appeared.

"Santa vaca," I said cleverly in Español, one of the native languages in
New Mexico, "you have proven my point! Look at the string of dry years
that ended in 1956."

"It sure looks that way, but how certain are you? Maybe a string of dry
years like that could happen any old time just by drawing randomly from
this population."

"Well," I said, "let's look at the distribution of values. Maybe that
will tell us something."

"Okay," he said, "but let's do some housekeeping first." He pointed and
clicked a few times, using the *File \| Save Workspace* command on R's
GUI menu to do something called saving the workspace. I thought it was
quite clever that the menu item was called *Save Workspace*, and
commented that the software might not be so shabby after all.

John said simply, "It's free."

Not one to concede a point easily, I said, "You sure have to type a lot.
I mean *pcp.df\$Annual* seems like a long way to refer to the annual
precip data."

"Oh, shut up," he said.

I could tell I was getting on his nerves. Maybe time to buy more
cerveza.

"Watch this," he said, and typed

attach(pcp.df) ; Annual

and there was a short listing that had alternating lines of years and
annual precip values. I could see that the values were listed right
under the associated years. Somehow R now knew that *Annual* referred to
*pcp.df\$Annual*.

"Of course," McBabble babbled, "that's what *attach* is for. You have to
be careful, though, that you don't attach something else that also has a
column named *Annual* in it or you will get confused." Forewarned is
forearmed, I was not about to do something so unequivocally irrational
(in Español this is called "loco").

By now, I had forgotten that we were going to look at the distribution
of annual precipitation values. McBabble, however, with a mind like a
steel sieve, typed

hist(Annual,breaks=2\*c(0:10))

and I saw a histogram that was a thing of beauty, nearly symmetrical and
as bell shaped as a cathedral bell. Noting my admiration, he added

x\<-0.2\*c(1:100)

curve(2\*86\*dnorm(x,mean=mean(Annual),sd=sd(Annual)),add=T,col=\"blue\")

and, although I now saw a normal distribution curve superimposed on the
histogram --- with a beautiful fit --- I sure didn't understand what he
had typed. I scratched my head.

He said, "It's very simple. I generated 100 values that range across the
*x*-axis values of the histogram, and then I used the mean and standard
deviation of *Annual* to calculate values from a normal distribution.
Then I added a blue (col="blue") curve on top of (add=T or add=TRUE) the
histogram."

"Why did you multiply by 2 \* 86?" I asked.

"The *dnorm* function generates probability densities," he said, looking
at me like I was the prime local manifestation of density, "so the area
under the whole curve of *dnorm* output is 1. If you look at our
histogram, the area is obviously the sum of the counts times the width
of the boxes. I multiplied by 5 and 86 to make the areas the same."

Obviously.

"Now," I said, "we've lost our time-series plot of the data. Can we look
at the time series together with this plot?

"R can do anything," he said.

"And it's free," I responded. I was catching on.

"Do you want the time series on top, or the histogram on top?" he asked.

"Histogram on top," I responded.

"Okay," he said, and typed

par(mfrow=c(2,1))

hist(Annual,breaks=2\*c(0:10),freq=T)

curve(2\*86\*dnorm(x,mean=mean(Annual),sd=sd(Annual)),add=T,col=\"blue\")

plot(Year,Annual,type="o")

"What's *par*?" I asked.

"It sets graphics parameters. This says plot in panels of two rows and
one column, going down rows before moving across columns."

This struck me as really weird notation, but I figured that if I could
look it up, it was probably all right.

"Can I find help on these things? They might be hard to remember."

McBabble said nothing. He just typed

help.start()

and his web browser started up. He was obviously in a hypertext help
system, which I knew I would need.

"Now look at those values in the decade before 1955. Don't they look low
to you?" I asked.

"Sure," he replied. "But don't you remember when you were an
*estudiante* (this is the New Mexican word for *victim*) in statistics
class? There you had to prove just how confident you were that something
in your data was anomalous."

Now that he was using Latin on me, I knew I was in trouble.

"First off," he continued, "you are talking about a ten-year average.
Let's figure out what the average was for that ten years and then think
about how it might compare to other ten-year averages that might be
drawn at random."

"What else could I possibly want to do this afternoon?" I said as I
signaled the waiter for another round. Lubrication would help.

"First of all, remember that some years are missing. Let's find out for
sure which rows in our data tables we are talking about." He typed

pcp.df\$Year,

peered at the screen for a moment, and then typed

pcp.df\$Year\[30:39\]

I looked at the screen and saw the numbers 1946 -- 1955 listed ---
exactly the years that we were interested in.

John typed

lowYrs\<-c(30:39)

lowYrAvg\<-mean(pcp.df\$Annual\[lowYrs\])

mean(Annual)

lowYrAvg

and then whistled.

I looked at the numbers: the average annual precip over the whole period
of record was 8.44 inches, and the average during our dry period was
6.12 inches. Looking at the histogram, I could see that 6.12 inches was
way up on the histogram shoulder --- actually near the top. So that did
not look very anomalous to me.

"Why do you whistle?" I asked.

"Because it's so anomalous," he answered.

"But 6.12 is near the center of the histogram," I replied.

"You forget that you are working with a ten year average, not just one
data point," he answered. "Because you are averaging 10 data points, you
have to squeeze the histogram by a factor of square root of 10, or 3.16"
he said. "When you do that in your mind, you can see that 6.12 is way
out on the tail of the histogram."

Cerveza fog was setting in.

"Why do you squeeze the histogram again?" I asked.

"Here, let me illustrate," he said. "Let's take 100 random samples, each
one with 10 values, out of this population."

"Why not," I responded. It was clear that I was not going anywhere else
that afternoon.

McBabble said, "I'll make an empty array of 100 rows of ten columns, and
then fill that with year subscripts (these are the row numbers in
pcp.df\$Year) selected at random." He typed

ranYr=array(dim=c(100,10));ranYr

for(iSamp in 1:100)ranYr\[iSamp,\]\<-sample(1:86,10);ranYr

"Notice," he said, "that the default for sample is to take a sample and
not put it back, so we don't have any years repeated in a row. Now let's
go from year subscripts to annual precip for the corresponding year."

ranAnnual\<-Annual\[ranYr\];ranAnnual

"Well, we had random year subscripts in 100 rows of 10 samples," he
noted, "but now we just have a thousand individual samples. Let's put
them back in rows of 10," and typed

dim(ranAnnual)\<-c(100,10);ranAnnual

"Notice how I can put new values into the dimension descriptor for an
array? I could have used c(10,10,10) and ended up with a 3-dimensional
array; the only thing I can't do is change the total size of the array.
Now let's take the means of these 10-year random samples and look at
their distribution."

tenYrMns\<-rowMeans(ranAnnual);tenYrMns

hist(tenYrMns)

"Now, what do you think of the value 6.12 for the mean of the years 1946
-- 1955?"

"Wow," I answered, "6.12 is way out on the tail of the distribution."

"Yes, it is," he said. "Let's compare the standard deviations of these
distributions."

sdAnnual\<-sd(Annual)

sd10YrMn\<-sd(tenYrMns)

sdAnnual;sd10yrMn;sdAnnual/sd10YrMn

"You see, the histogram curve of 10-year means is the histogram curve of
single years, squeezed by the square root of 10, which is 3.16. Of
course our amount of squeeze in our sample is a little bit different
from this because random samples vary a little bit from the ideal, but
you get the idea that the distribution of sample means is always
narrowed by the square root of the sample size."

"How can I remember that when I'm adjusting standard deviations for
means of *n* samples I should divide by square root of *n* and not
(*n*-1)?" I asked, remembering that in the formula for calculating
standard deviations I divide the total of the squared deviations by
(*n*-1).

"Easy," he responded. "Just think of a single sample as an average with
an *n* of 1. You can't very well divide by the square root of zero, so
the formula has to use *n*, not (*n*-1)."

"Yeah, easy," I responded. "For you to say," I thought.

"So do you always have to go through this generating random samples and
looking at the distribution of their means?" I asked.

"Oh, no," he replied. "I just did that to prove the point that the
distribution gets squeezed. Once you trust that that is how it works,
you just use the formula for *z*. What *z* does is tell you directly how
many standard deviations you are away from the mean, where the standard
deviation is adjusted for the number of samples that you have averaged
together. Here's how you do it."

z55\<-(mean(Annual\[lowYrs\])-mean(Annual))/(sd(Annual)/sqrt(10));z55

"And this is the number of standard deviations that the average annual
precip for the 10-year period ending in 1955 was below the long-term
average annual precip."

"Wow, so it was 2.7 standard deviations below normal," I said.

"That's right, so it was pretty abnormal. But how much confidence do you
have in its abnormality?"

I wasn't sure I knew what he meant. Suspecting that he was going to make
me calculate my confidence, I asked, "How can we figure that out?"

"Here's what we'll do," he said. "We'll divide the probability
distribution into two parts. I want you to assert that the result was
abnormal --- not the result of chance --- and then we'll figure out what
part of the distribution you own. Okay?"

"Why not?" I answered, not having a clue what was going to happen.

"Okay, so you say that your result, 2.7 standard deviations
[below]{.underline} the mean, was abnormal. What would you say about a
result that was 2.7 standard deviations [above]{.underline} the mean?"

"Abnormal," I said.

"Okay, and what about something that was 3 standard deviations below the
mean?"

"Abnormal," I said.

"You can probably see where I'm going with this," said McBabble. If you
say 2.7 standard deviations below the mean was abnormal, then you have
to admit that everything in both tails of the distribution beyond 2.7
standard deviations is also abnormal."

"Yeah, I guess so. How much of the distribution is that?" I asked.

McBabble typed

pnorm(-2.68,mean=0,sd=1)\*2

and said, "Zero point 7 percent."

"Okay," I responded, "I own 0.7% of the distribution. What does that
mean?"

"It means that all the things that you define as abnormal --- and this
has to include the tails of the distribution that have values even more
abnormal than the one that you have in your sample --- might have
happened by chance 0.7% of the time. By saying your value is abnormal,
you are doing what we call *rejecting the null hypothesis*. In
statistics, the null hypothesis is always that the result happened by
random chance."

"So you're saying that my definition of *abnormal* could have happened
by chance only 0.7% of the time with randomly selected data. That's not
very much."

"That's right," said McBabble. "You can be pretty confident that your
assertion is right that the decade 1946 -- 1955 was pretty weird. In
statistics, we call the distribution tail area that we calculated the
*p* value. It is the chance that you are wrong in saying that your value
is an anomaly, or saying that you reject the null hypothesis. Because
your *p* value is less than 1%, you are more than 99% confident that
your definition of *abnormal* is not the result of chance."

"I knew it was a weird decade," I said.

"You did not. You had a qualitative hunch. And you still don't know it,
but now you know what the chances are that you are wrong."

"One other thing. We did all our analysis assuming that annual precip
values are drawn at random. Isn't it likely that one dry year is
followed by another, so the probabilities change if we know the prior
year's precip?"

"You are absolutely right," he responded. "But that gets into temporal
correlation, and that's a topic for another discussion, and my batteries
R getting low."

"Well I sure was reminded that statistics R fun," I responded.

"Statistics R free," he said as he shut his laptop and stood to go.

"Right," I thought as I picked up the tab.

Raymond D. Watts

USGS Fort Collins Science Center and Rocky Mountain Geographic Science
Center