Posts filed under 'Good Statistics'
Here is the link.
I’ve decided to jump ahead a few chapters. Chapters 10 - 13 are very important and cover material that comprises 90% of all the actual statistics that is practiced by civilians. Topics like “testing” and regression—how they are done in classical and Bayesian statistics, why these methods are too sure of themselves, and why observable statistics is the only proper way.
But I can tell that I am testing the patience of my audience, so I will leave these more technical chapters for the book itself.
Thus, here I return to something eminently practical: HOW TO CHEAT WITH STATISTICS.
It is important these days for people to know how to get away with as much as they possibly can. This chapter shows you how to do it.
There are no cheap methods like data fudging or just plain lying—those techniques are for pikers. No: what I give you is genuine, sophisticated gold. Tricks you can actually use and get away with. Tricks that work.
I must be out of my mind to give these secrets away for free, but it is a measure of how much I love you, my audience, my faithful readers.
Only an excerpt is in this posting. To get the whole Chapter, you’ll have to download it. Here is the link.
2. Conditioning
A typical academic study is one, say, that gathers two groups of college kids, maybe about 50 in each set, and has them do some task or asks them to rate something. Another study gathers data from a small area, say a neighborhood in a city, where the sample size may be as high as a few hundred, and asks sociological and economic questions of the people that live there. A medical study might try two treatments in two groups of a hundred or so people. When the data from these studies are in, the results are compiled and papers are published. Claims are made in these papers. The college kids paper will say that people act one way and not another; the city paper will say that poor people have less money; and the medical paper will claim treatment A is better than treatment B.
We already know that if all these researchers wanted to do was to say something about their datasets, then they do not need statistics or probability models. They can look at their data and say, yes, more people got better under treatment A than under treatment B. They would be finished. Evidently, the creators of these studies do not want to make statements only about the past data; they want to imply their findings are more widely applicable.
By far the majority of these kinds of studies, published in academic journals, concern humans. As of this writing, there are over 6.6 billion humans alive, about 100 billion are dead, and God only knows how many more are yet to live. Incidentally, whatever you do, do not mention these facts in your results (unless, of course, you happen to be writing about demography), it will weaken your argument.
Are the results from the college kids study applicable to all humans? All those that lived in the past, those that will live in the future, even those that live now but not in the town in which the college lie? Those who are in their 50s?, 80s? who are less than 10? Poorer people and those with enough money to “get a degree”? (Kids go to college to “get a degree” nowadays, and not usually for anything else. Well, maybe socialization. These are rational choices given the way things are.) Kids at other universities? Let’s be clear: the researchers will gather data on their 100 kids, create a probability model, and since they have read this book, they will not just make a statement about the parameters, but calculate the probability distribution of future observables. The only problem is, about whom do we apply this probability distribution?
Before we answer that, think about the medical trial, which was conducted at a hospital in a city on the East Coast of the United States of America. The physicians also use their data to create a probability distribution of future patients. But who exactly are these patients? People who live in other cities on the east coast?, anywhere in the USA? Canada, too? Or only cities of a certain size? Or do the future patients merely have to “look like” the patients in the old data; that is, be of the same ages, sex ratio, weights, economic condition, have eaten the same things in their lifetimes, traveled to the same places, engaged in the same activities, and so on. Would it have applied to the people who used to be alive, and to people not yet born, indefinitely into the future?
Nobody knows the answers to these questions, which is highly in your favor, especially if you have just completed a study using data “at hand”, that is, that was easy for you to collect. You certainly want to imply that your results are as broadly applicable as possible because this makes you more of an expert than somebody who merely claims to know the habits of a small group of college kids in the year 2008 only, in city C and who are unmarried, between 19 and 22 years old, and whose parents are upper middle class, etc. Openly stressing these limitations might be noble and correct, but it will not get you far. State your results in terms of all people. For example, say “People choose option A over B which gives weight to our theory of psychology.” Do not say, “College kids in our freshman psychology class, who might not be anything like the rest of the population, carried out an experiment for us—and surely they took this task seriously—and…”
Same thing in the medical trial. Emphasize your small p-value, spend more time talking about how the two groups of patients (those that received treatment A and those that got B) were not different than one another. Tell how there were roughly equal numbers of men and women in both treatments, and the same with age, weight, etc. This is an excellent strategy because it is useful information: if the two groups did differ, then your results may be biased. Well, this is a wonderful distraction because it allows you to ignore or downplay the discussion of how your results might only be useful for a small subset of patients.
In short, be sloppy describing the nature of your sample; or, rather, say as much about your sample as you like, but say little or nothing about whom you expect your results are applicable. Certainly imply that all humanity falls under your results. With any luck, a reporter will find your paper and help you along this road by summarizing your results, leaving out all hint of limitation with his headline, “Kumquats reduce risk of toenail cancer.”
5. Publishable p-values
Most journals, say in medicine or those serving fields ending with “ology”, are slaves to p-values. Papers have a difficult, if not impossible, time getting published unless authors can demonstrate for their study a p-value that is publishable, that is, that is less than 0.05. Sometimes, the data are not cooperative and the p-value that you get from using a common statistic is too large to see the light of print. This is bad news, because if you are an academic, you must publish papers else you can’t get grants, and if you don’t get grants, then you do not bring money into your university, and if you don’t bring money into your university, you do not get tenure, and if you do not get tenure, then you are out the door and you feel shame.
So small p-values are important. I of course advise against using classical statistics methods, but if you are forced to (and some journal editors insist on it), then all is not lost if an initial large p-value is found. In fact, I would go so far to say that if you cannot find a publishable p-value in any situation, then you are not trying hard enough. There are several ways to lower your p-value.
The most well known is to increase your sample size. This one is a lock. Let’s take a look at the t-test statistic from Chapter 10 to see why.
(see the book)
There is a mathematical phrase that begins “without loss of generality” which I now invoke by letting, for ease of notation, nA = nB = n and s2 = s2 = s2 , so that t(x) becomes
(see the book)
Remember that we want a large statistic, a large t, the larger the better, because larger ts mean smaller p-values. Do you see the trick? A larger n means a larger t! All you have to do is to increase your sample size and just wait for the small p-values to start rolling in. This trick always works in any classical situation, even when the difference xA − xB is too small to be of interest to anybody. This is why having a small p-value is called attaining statistical significance and not practical or useful significance.
Incidentally, this trick also works in Bayesian statistics in the sense that the posterior distribution of μ A − μ B will have most probability above or below zero. But it fails miserably in modern observable statistics because a trivial difference in μ A − μ B won’t make a tinker’s dam worth of difference in the probability distribution of future observables.
The next trick, if you cannot increase your sample size, is to change your statistic. This comes from the useful loophole in classical theory that there is no rule which specifies which statistic you must use in any situation. Thus, though some creativity and willingness to spend time with your statistical software, you can create small p-values where other people see only despair. This isn’t so easy to do in R because you have to know the names of the alternate statistics, but it’s cake in SAS, which usually prints out dozens of statistics in standard cases, which is one reason SAS is worth its exorbitant price. Look around at the advertising brochures of statistical software and you will see that the openly boast of the large number of tests on offer.
For example, for use in “testing differences between proportions”, just off the top of my head I can think of the z statistic, the proportions test with and without correction for continuity (two or three to choose from here), chi-squared test, Fisher’s exact test, McNemar’s test, logistic regression. There are dozens more and teams of academic statisticians constantly add to the pile. Don’t believe it? Here’s a small table of these tests for the TSD/Sex data from Chapter 11.
| Test |
p-value |
| Prop test |
0.78 |
| Fisher’s |
0.70 |
| Logistic Reg. |
0.52 |
| chi-squared |
0.50 |
| z test |
0.49 |
| McNemar’s |
0.24 |
Because I was only able to get to 0.24 just means I didn’t try hard enough. Which is the correct p-value? They all are; that’s the beauty of this trick. Not one of these p-values is more “right” than any other one. Each is valid. If all you know is classical statistics, let this knowledge sink in. It should prove to you that p-values are not what you probably thought they were.
For ‘testing differences between means”, there is the t-test (a couple of versions of this, actually), Wilcox test (also called Mann- Whitney), sign tests, Spearman correlation tests, Kendall’s tau, Kruskal-Wallis test, Kolmogorov-Smirnov test, permutation test, Friedman two-way analysis of variance—I’m running out of breath—and many more. Here’s some of those tests for the advertising data:
| Test |
p-value |
| Spearman |
0.87 |
| Perm. |
0.20 |
| t-test |
0.19 |
| Wilcox |
0.14 |
| Kol.-Smi. |
0.08 |
Nearly there!
Please remember that in this example, like the previous one, the data is the same; the only thing that changes is that classical statistical test.
The key to this deceit is to never admit what you did. When it comes time to write up your result boldly and authoritatively state, “We used Johnston’s (Johnston, 1983) frammilax test for differences in means.” Tossing in a citation always cows potential critics; tossing in two or more guarantees editorial acquiesence. Do not tell the reader that you went through a dozen tests to find the lowest p-value. Act as if “Johnston’s test” was what you had in mind all along.
This technique is unavailable in Bayesian or observable statistics. True, you can change your default prior distribution on the parameters or even change the model (see below), but editors in most fields are still suspicious of modern methods and tend to be conservative and will likely insist on a well-known default. There will be more room for creativity in, say, ten years when modern methods become familiar.
Our last option, if you cannot lower your p-value any other way, is to change what is accepted as publishable. So, instead of a p-value of 0.05, use 0.10 and just state that this is the level you consider as statistically significant. I haven’t seen any other number besides 0.10, however, so if your p-value is larger than this the best you can do is to claim that your results are “suggestive” or “in the expected direction.” Don’t scoff, because this sometimes works.
You can really only get away with this in secondary and tertiary journals (which luckily are increasing in number) or in certain fields where the standard of evidence is low, or when your finding is one which people want to be true. This worked for second-hand smoking studies, for example, and currently works for anything negatively associated with global warming.
June 5th, 2008
Here is the link.
The going started getting tough last Chapter. It doesn’t get any easier here. But stick with it, because once you finish with this Chapter you will see the difference between classical/Bayesian and modern statistics.
Here is the gist:
- Start with quantifying your uncertainty in an observable using a probability distribution
- The distribution will have parameters which you do not know
- Quantify your uncertainty in the parameters using probability
- Collect observable data ,which will give you your updated information about the parameters which you still do not know and which still have to be quantified by a probability distribution
- Since you do not care about the parameters, and you do care about future observables, you quantify your uncertainty in these future observables given the uncertainty you still have in the parameters (through the information present in the old data).
If you stop at the parameters, step 4, then you are a regular Bayesian, and you will be too certain of yourself.
This Chapter shows you why. The computer code mentioned in the homework won’t be on-line for a week or so. Again, some of you won’t be able to see all Greek characters, and none of the pictures are given. You have to download the chapter. Here is the link.
CHAPTER 9
Estimating and Observables
1. Binomial estimation
In the 2007-2008 season, the Central Michigan football team won 7 out of 12 regular season games. How many games will they win in the 2008-2009 season? In Chapter 4, we learned to quantify the probability in this number using a binomial distribution, but we assumed we knew p, the probability of winning any single game. If we do not know p, we can use the old data from last season to help us make a guess about its value. It helps to think of this old data as a string of wins and losses. So that, for the old x, we saw x1 = 0, x2 = 1, . . . , x12 = 1, which we can summarize by k = i xi , where k = 7 is the total number of wins in n = 12 games.
Here’s the binomial distribution written with an unknown parameter
(see the book)
where θ is the success parameter and k the number of successes we observed out of n chances.
How do we estimate θ? Two ways again, a classical and a modern. The classical consists of picking some function of the observed data and calling it θ, and then forming a confidence interval. In R we can get both at once with this function
binom.test(7,12)
where you will see, among other things (ignore those other things for now),
95 percent confidence interval:
0.2766697 0.8483478
sample estimates:
probability of success
0.5833333
This means that θ = 0.58 = 7/12 so again, the estimate is just the arithmetic mean. The 95% confidence interval is 0.28 to 0.84. Easy. This confidence interval has the same interpretation as the one for the μ, which means you cannot say there is a 95% chance that θ is in this interval. You can only say, “either θ is in this interval or it is not.”
Here is Bayes’s theorem again, written as functions like we did for the normal distribution
(see the book)
We know p(k|n, θ, EB ) (this is the binomial distribution), but we need to specify p(θ|EB ), which describes what we know about the success parameter before we see any data, given only EB (p(k|n, EB ) will pop out using the same mathematics that gave us p(x|EN ) in equation (17)). We know that θ can be any number between 0 and 1: we also know that it cannot be exactly 0 or 1 (see the homework). Since it can be any number between 0 and 1, and we have no a priori knowledge which number is more likely than any other, it may be best to suppose that each possible value is equally likely. This is the flat prior again (1Like before, there are more choices for this prior distribution, but given even a modest sample size, the differences in the distribution of future observables due to them is negligible). Again, technically EB should be modified to contain this information. After we take the data, we can plot p(θ|k, n, EB ) and see the entire uncertainty in θ, or we can pick a “best” value, which is (roughly) θ = 0.58 = 7/12, or we can say that there is a 95% chance that θ is in the (approximate) interval 0.28 to 0.84. I say “roughly” and “approximate” here, because the classical approximation to the exact Bayesian solution isn’t wonderful for the binomial distribution when the sample size is small. The homework will show you how to compute the precise answers using R.
2. Back to observables
In our hot little hands, we now have an estimate of θ which equals about 0.58. Does this answer the question we started with?
That question was How many games will CMU win in the 2008-2009 season? Knowing that θ equals something like 0.58 does not answer this. Knowing that there is a 95% chance that θ is some number between 0.28 to 0.84 also does not answer the question. This question is not about the unobservable parameter θ, but about the future (in the sense of not yet seen) observable data. Now what? This is one of the key sections in this entire book, so take a steady pace here.
Suppose θ was exactly equal to 0.58. Then how many games will CMU win? We obviously don’t know the exact number even if we knew θ, but we could calculate the probability of winning 0 through 12 games using the binomial distribution, just as we did in Chapters 3 and 4. We could even draw the picture of the entire probability distribution given that θ was exactly equal to 0.58. But θ might not be 0.58, right? There is some uncertainty in its value, which is quantified by p(θ|kold , nold , EB ), where now I have put the subscript “old” on the old data values to make it explicit that we are talking about the uncertainty in θ given previously observed data. The parameter might equal, say, 0.08, and it also might equal 0.98, or any other value between 0 and 1. In each of these cases, given that θ exactly equalled these numbers, we could draw a probability distribution for future games won, or knew given nnew = 12 (12 games next season) and given the value of θ.
Let us draw the probability distribution expressing our uncertainty in knew given nnew = 12 (and EB ) for three different possible values of θ.
(see the book)
If θ does equal 0.08, we can see that the most likely number of games next season is 1. But if θ equals 0.58, the most likely number of games won is 7; while if θ equals 0.98, then CMU will most likely win all their games.
This means that the picture on the far left describes our uncertainty in knew if θ = 0.08. What is the probability that θ = 0.08? We can get it from equation (19), from p(θ|kold = 7, nold =12, EB ). The chance of θ = 0.08 is about 1 in 100 million (we’ll learn how the computer does these calculations in the homework). Not very big! This means that we are very very unlikely to have our uncertainty quantified by the picture on the left. What is the chance that θ = 0.98? About 3 in a trillion! Even less likely. How about 0.58? About 3 in 10,000. Still not too likely, but far more likely than either of those other values.
We could go through the same exercise for all the other values that θ could take, each time drawing a picture of the probability distribution of knew . Each one of these would have a certain probability of being the correct probability distribution for the future data, given that its value of θ was the correct value. But since we don’t know the actual value of θ, but we do know the chance that θ takes any value, we can take a weighted sum of these individual probability distributions to produce one overall probability distribution that completely specifies our uncertainty in knew given all the possible values of θ. This will leave us with
(see the book)
Stare at equation (20) for two minutes without blinking. This, in words, is the probability distribution that tells us everything we need to know about future observables knew given that we know there will be nnew chances for success this year, also given that we have seen the past observables kold and nold , and assuming EB is true. Think about this. You do not know what future values of k will be, do you? You do know what the past values are, right? So this is the way to describe your uncertainty in what you do not know given what you do know, taking full account of the uncertainty in θ, which is not of real interest anyway.
The way to get to this equation uses math that is beyond what we can do in this class, but that is unimportant, because the software can handle it for you. This picture shows you what happens. The solid lines are the probability distribution in equation (20). The circles plotted over it are the probability distribution of a regular binomial assuming θ exactly equals 0.58. The key thing to notice is that the circles distribution, which assumes θ ≡ 0.58 is too tight, too certain. It says the center values of 6 to 8 are more certain than is warranted (their probability is higher than the actual distribution). It agrees, coincidentally only, with the probability that the future number of wins will be 5 or 9, but then gives too little probability for wins less than 5 or greater than 9.
The actual distribution of future observable data (20) will always be wider, more diffuse and spread out, less certain, than any distribution with a fixed θ. This means we must account for uncertainty in the parameter. If we do not, we will be too certain. And if all we do is focus on the parameter, using classical or Bayesian estimates, and we do not think about the future observables, we will be far, far more certain than we should be.
3. Even more observables
Let’s return to the petanque example and see if we can do the same thing for the normal distribution that we just did for the binomial. The classical guess of the central parameter was μ = −1.8 cm, which matches the best guess Bayesian estimate. The confidence/credible interval was -6.8 cm to 2.8 cm. In modern statistics, we can say that there is a 95% chance that μ is in this interval. We also have a guess for σ, and a corresponding interval, but I didn’t show it; the software will calculate it. We do have to think about σ as well as μ, however—both parameters are necessary to fully specify the normal distribution.
As in the binomial example, we do not know what the exact value of (μ, σ) is. But we have the posterior probability distribution p(μ, σ|xold , EN ) to help us make a guess. For every particular possible value of (μ, σ), we can draw a picture of the probability distribution for future x given that that particular value is the exact value.
(see the book)
The picture shows the probability densities for xnew for three possible values of (μ, σ). If (μ = −6.8 cm, σ = 4.4 cm), the most likely values of xnew are around 10 cm, with most probability given to values from -20 cm to 0 cm. On the other hand, if (μ = 2.8 cm, σ = 8.4 cm), the most likely values of new x are a little larger than 0 cm, but with most probability for values between -20 cm and 30 cm. If (μ = −1.8 cm, σ = 6.4 cm), future values of x are intermediate of the other two guesses. These three pictures were drawn (using the Advanced code from Chapter 5) assuming that the values of (μ, σ) are the correct ones. Of course, they might be the right values, but we do not know that. Instead, each of these three guesses, and every other possible combination of (μ, σ), has a certain probability, given xold , of being true.
Given the old data, we can calculate the probability that (μ, σ) equals each of these guesses (and equals every other possible combination of values). We can then weight each of the new x distributions according to these probabilities and draw a picture of the distributions of new values given old ones (and the evidence EN ) like we just did for the binomial distribution. This is
(see the book)
Here is a picture of this distribution (generated by the computer, of course)
(see the book)
The solid line is equation (21), and dashed is a normal distribution with (μ = −1.8 cm, σ = 6.4 cm). The two distributions do not look very different, but they certainly are, especially for very large or very small values of xnew . The dashed line is too narrow, giving too much probability for too narrow a range of xnew . In fact, for distribution (21), values greater than 10 cm are from the true distribution are twice as likely as the normal distribution where we plugged in a single guess of (μ, σ); values greater than 20 cm are six times as likely. The same thing is repeated for values less than -10 cm, or less than -20 cm, and so on. Go back and read Chapter 6 to refamiliarize yourself with the fact that very small changes in the central or variance parameter can cause large changes in the probability of extreme numbers.
The point again, like in the binomial example, is that using the plug-in normal distribution, the one where you assume you know the exact value of (μ, σ), leads you to be far more certain than you really should be. You need to take full account of the uncertainty in your guesses of (μ, σ), only then will you be able to full quantify the uncertainty in the future values xnew .
May 30th, 2008
Here is the link.
This is where it starts to get complicated, this is where old school statistics and new school start diverging. And I don’t even start the new new school.
Parameters are defined and then heavily deemphasized. Nearly all of old and new school statistics entire purpose is devoted to unobservable parameters. This is very unfortunate, because people go away from a parameter analysis far, far too certain about what is of real interest. Which is to say, observable data. New new school statistics acknowledges this, but not until Chap 9.
Confidence intervals are introduced and fully disparaged. Few people can remember that a confidence interval has no meaning; which is a polite way of saying they are meaningless. In finite samples of data, that is, which are the only samples I know about. The key bit of fun is summarized. You can only make one statement about your confidence interval, i.e. the interval you created using your observed data, and it is this: this interval either contains the true value of the parameter or it does not. Isn’t that exciting?
Some, or all, of the Greek letter below might not show up on your screen. Sorry about that. I haven’t the time to make the blog posting look as pretty as the PDF file. Consider this, as always, a teaser.
For more fun, read the chapter: Here is the link.
CHAPTER 8
Estimating
1. Background
Let’s go back to the petanque example, where we wanted to quantify our uncertainty in the distance x the boule landed from the cochonette. We approximated this using a normal distribution with parameters m = 0 cm and s = 10 cm. With these parameters in hand, we could easily quantify uncertainty in questions like X = “The boule will land at least 17 cm away” with the formula Pr(X|m = 0 cm, s = 10 cm, EN ) = Pr(x > 17 cm|m = 0 cm, s = 10 cm, EN ). R even gave us the number with 1-pnorm(17,0,10) (about 4.5%). But where did the values of m = 0 cm and s = 10 cm come from?
I made them up.
It was easy to compute the probability of statements like X when we knew the probability distribution quantifying its uncertainty and the value of that distribution’s parameters. In the petanque example, this meant knowing that EN was true and also knowing the values of m and s. Here, knowing means just what it says: knowing for certain. But most of the time we do not know EN is true, nor do we know the values of m and s. In this Chapter, we will assume we do in fact know EN is true. We won’t question that assumption until a few Chapters down the road. But, even given EN is true, we still have to discern the values of its parameters somehow.
So how do we learn what these values are? There are some situations where are able to deduce either some or all of the parameter’s values, but these situations are shockingly few in number. Nearly all the time, we are forced to guess. Now, if we do guess—and there is nothing wrong with guessing when you do not know—it should be clear that we will not be certain that the values we guessed are the correct ones. That is to say, we will be uncertain, and when we are uncertain what do we do? We quantify our uncertainty using probability.
At least, that is what we do nowadays. But then-a-days, people did not quantify their uncertainty in the guesses they made. They just made the guesses, said some odd things, and then stopped. We will not stop. We will quantify our uncertainty in the parameters and then go back to what is of main interest, questions like what is the probability that X is true? X is called an observable, in the sense that it is a statement about an observable number x, in this case an actual, measurable distance. We do not care about the parameter values per se. We need to make a guess at them, yes, otherwise we could not get the probability of X. But the fact that a parameter has a particular value is usually not of great interest.
It isn’t of tremendous interest nowadays, but again, then-a-days, it was the only interest. Like I said, people developed a method to guess the parameter values, made the guess, then stopped. This has led people to be far too certain of themselves, because it’s easy to get confused about the values of the parameters and the values of the observables. And when I tell you that then-a-days was only as far away as yesterday, you might start to be concerned.
Nearly all of classical statistics, and most of Bayesian statistics is concerned with parameters. The advantage the latter method has over the former, is that Bayesian statistics acknowledges the uncertainty in the parameters guesses and quantifies that uncertainty using probability. Classical statistics—still the dominate method in use by non-statisticians1—makes some bizarre statements in order to avoid directly mentioning uncertainty. Since classical statistics is ubiquitous, you will have to learn these methods so you can understand the claims people (attempt to) make.
So we start with making guesses about parameters in both the old and new ways. After we finish with that, we will return to reality and talk about observables.
2. Parameters and Observables
Here is the situation: you have never heard of petanque before and do not know a boule from a bowl from a hole in the ground. You know that you have to quantify x, which is some kind of distance. You are assuming that EN is true, and so you know you have to specify m and s before you can make a guess about any value of x.
Before we get too far, let’s set up the problem. When we know the values of the parameters, like we have so far, we write them in Latin letters, like m and s for the Normal, or p for the binomial. We always write unknown and unobservable parameters as Greek letters, usually μ and σ for the normal and θ for the binomial. Here is the normal distribution (density function) written with unknown parameters:
(see the book)
where μ is the central parameter, and σ 2 is the variance parameter, and where the equation is written as a function of the two unknowns, N(μ, σ). This emphasizes that we have a different uncertainty in x for every possible value of μ and σ (it makes no difference if we talk of σ or σ 2 , one is just the square root of the other).
You may have wondered what was meant by that phrase “unobservable parameters” last paragraph (if not, you should have wondered). Here is a key fact that you must always remember: not you, not me, not anybody, can ever measure the value of a parameter (of a probability distribution). They simply cannot be seen. We cannot even see the parameters when we know their values. Parameters do not exist in nature as physical, measurable entities. If you like, you can think of them as guides for helping us understand the uncertainty of observables. We can, for example, observe the distance the boule lands from the cochonette. We cannot, however, observe the m even if we know its value, and we cannot observe μ either. Observables, the reason for creating the probability distributions in the first place, must always be of primary interest for this reason.
So how do we learn about the parameters if we cannot observe them? Usually, we have some past data, past values of x, that we can use to tell us something about that distribution’s parameters. The information we gather about the parameters then tell us something about data we have not yet seen, which is usually future data. For example, suppose we have gathered the results of hundreds, say 200, of past throws of boules. What can we say about this past data? We can calculate the arithmetic mean of it, the median, the various quantiles and so on. We can say this many throws were greater than 20 cm, this many less. We can calculate any function of the observed data we want (means and medians etc. are just functions of the data), and we can make all these calculations never knowing, or even needing to know, what the parameter values are. Let me be clear: we can make just about any statement we want about the past observed data and we never need to know the parameter values! What possible good are they if all we wanted to know was about the past data?
There is only one reason to learn anything about the parameters. This is to make statements about future data (or to make statements about data that we have not yet seen, though that data may be old; we just haven’t seen it yet; say archaeological data; all that matters is that the data is unknown to you; and what does “unknown” mean?). That is it. Take your time to understand this. We have, in hand, a collection of data xold , and we know we can compute any function (mean etc.) we want of it, but we know we will, at some time, see new data xnew (data we have not yet seen), and we want to now say something about this xnew . We want to quantify our uncertainty in xnew , and to do that we need a probability distribution, and a probability distribution needs parameters.
The main point again: we use old data to make statements about data we have not yet seen.
3. Classical guess
We need to find some way to map our evidence E and the past values of x into information about the parameters. There are lots of different ways to guess at parameter values, some easy and some hard, and these all fall into two broad classifications: yes, a classical and a modern.
We have past values of x and we want to know about future, or at least other, unknown values of x. Our evidence is E, which at least means that we know the probability distribution (Normal, say) of the observables. In this book we will also assume that E also means that knowledge of each individual observation is irrelevant to knowing what each other observation with be. We have to find a way to guess, or estimate, these unknown and unobservable parameters given E and the old data xold .
The classical way to do this is to pick an ad hoc function of the old data and label it f (xold ) = μ, where that “hat” indicates that the value of μ is only a guess. Most classical estimates have the goal that the estimate is “unbiased”, or Ex (μ − μ) = Ex (μ − f (xold )) = 0, meaning that the expected distance between the actual value of μ and the guess μ is 0. Sounds like a nice thing to have, unbiasedness, and it surely isn’t a bad idea, but it turns out to cause a lot of problems, most of which I cannot tell you about without introducing a lot of math. However, this criterion is not compelling because of that expected value business. Expected value with respect to what? Well, with respect to an infinite number of future (not yet observed) data x…which is just the data that we are trying to quantify the uncertainty of. Anyway, in R, to estimate the parameters of a normal distribution classically is easy, and you already know how to do it! If x is our old, previously observed data, x1 , x2 , …, xn , then
μ = mean(x) σ = sd(x)
The mean you already know to calculate. It is often written bar(x), and called “x bar”. When you see a data value with a bar over it, you know it is a mean. The observed variance of old data is (see the book), and the observed standard deviation of old data is the square root of that. Look at the formula and notice that the standard deviation is a measure of how far, on average, the old data values are away from the observed mean. The square is taken, (xi − bar(x))^22 , so that data values that were lower than the observed mean are treated the same as data values that were higher. (If you have missing data in x, recall Chapter 7, where we had to modify the function like this mean(x, na.rm=T; same for the sd function).
We’ll never calculate the observed standard deviation by hand. But it’s pretty convenient to have the observed mean stand in for our guess of μ. Unfortunately, because μ = mean, a lot of people have taken to calling μ (without the hat) the mean, which it most assuredly is not. μ is an unobservable parameter, while the mean is just the weighted sum of a bunch of data we have already observed. This is a subject that I’ll return to later.
4. Confidence intervals
OK, it might have been hard to understand all this so far, but it’s about to get weird, so be steady. The value μ we got before was precise; it is a known, observed number (it is the mean). But do we really believe, given the data and other evidence, that the exact, all-time, incorruptible, immutable value of μ is, to as many decimal places as you like, equal to μ? You may have guessed, by the subtle way I’ve asked that question, that the answer is “no.” And you’d be right! Suppose μ = 5.41. Maybe μ is 5.41, but it might also be, say, 5.40, or 5.39, or other values close by, mightn’t it? This is a fancy way to state that we are uncertain what the value of μ is. How do we express this uncertainty? Use probability? No. It is forbidden to use probability to quantify the uncertainty of parameter values in classical statistics.
Instead, classical statisticians use something called a confidence interval, which is an interval on the order of μ ± c(n), where c(n) is some number that usually depends on the number n of your data points and on the old data itself. Bigger c(n) lead to wider intervals; smaller c(n) lead to narrower ones. So you might expect that when you say that “I think μ is 5.41 plus or minus 4″ you have a better chance of being right then when you say “I think μ is 5.41 plus or minus 1″, because the former interval allows you greater scope of covering the actual (unobservable) value of μ. And, classically, you’d be dead wrong.
Which is why confidence intervals are one of the screwiest things to come out of the classical tradition, in that they fail utterly to do what they set out to do. But their use is so ubiquitous not to say iniquitous) that I’m afraid you are going to have to learn to interpret them. And they are one of the most important things you must learn in this book! because you will see confidence intervals everywhere, thus it is imperative you learn what they are and what they are not.
Part of the problem is that you simply cannot learn what a confidence interval is by reading most introductory statistics books. Take, for example, the very typical book Statistics: Informed Decisions Using Data by Sullivan (2007, pp. 448-449), often used in Stats 101 courses. He officially defines a confidence interval for an unknown parameter as “an interval of numbers” (p. 449), which is as pure a tautology as you’re ever likely to meet, and being a tautology, it is therefore, of course, true, but of no help (it says the confidence interval is an interval). But a page earlier, we find Sullivan implying that smaller intervals give us less confidence in the value of the parameter than larger intervals. This implication is, as I said above, false, and is no part of the actual, mathematical definition of a confidence interval.
Maybe something like this is more accurate:
[A] 95% level of confidence…implies that, if 100 different confidence intervals are constructed…we will expect 95 of the intervals to include the parameter and 5 to not include the parameter [p. 449].
Actually, we can expect nothing like this. And though this definition is closer to the truth, it is still false (to find out why, keep reading). Incidentally, classical theory lets you calculate confidence intervals at any level you want, but the only one you ever really see is the 95% interval, so that one is all I will talk about.
Here’s the actual definition. Suppose you gather some data and construct a confidence interval using the formula C1 = {μ ± c(n)} (the actual formula is not of much interest to us; the software will give us the interval automatically). That is, C1 is the interval calculated using the data we collected. Now imagine (incidentally, this is all you can do) that you re-collect your data in exactly the same way, where every physical thing is exactly the same as it was when you collected it the first time. That is, the state of the universe has to be identical to where it was when you first collected your data. Except that it must be “randomly” different, or different in ways that you know nothing about. Very well, you now have a second data set equal in every way to the first, except that it is “randomly” different, whatever that means. You then construct a new confidence interval C2 using the exact same formula on this second set of data (which is also the same size, n). Now do it all again and construct C3 , and again for C4, and again and again an infinite number of times. When you are done, 95% of those intervals will cover the actual value of μ.
(see the book)
This is shown in the picture for the first eight confidence intervals (this is all simulated data). The true value of μ is indicated by the solid line. Some of the intervals “cover”, i.e. contain, the true value of μ, and some do not. More than that, we cannot say. Our confidence interval, the bottom bold one, is the only confidence interval we’ll actually see; the others are hypothesized entities that are conjured into existence if confidence intervals are properly interpreted.
I only showed the first 8 (out of an infinite number of) confidence intervals (that must exist for every problem you ever do). If you only repeat your experiment a finite number of times, and therefore only have a finite number of confidence intervals, say, 1,000,000, then it is false that we expect any number of them will cover the true value of μ: stopping constructing confidence intervals at any finite value invalidates the interpretation that 95% of intervals will cover the actual value of μ.
Yes, this is the actual definition, but saying it this way leaves a bad taste in one’;s mouth, especially because of that bit about “infinite” numbers of repetitions. Statisticians, feeling uneasy about infinities, and their physical impossibility, usually resort to the euphemism “long run” to describe the number of repetitions needed. They know very well that, mathematically, long run equals infinite, but saying “in the long run” gives the comfortable impression that all you need is a lot, and not an infinite number, of repetitions.
By now you are thinking, “OK, I get it. So what? What you’re saying is just a quibble. Who cares about infinities or long runs, anyway. Give me some information I can use! What do you do with your confidence interval, the one you just constructed? What does it mean?”
Nothing. Not a thing. It certainly does not mean that you are 95% sure that your interval contains the actual value of μ. That is, you cannot, under any circumstances, say that “There is a 95% chance that the true value of μ lies in the 95% confidence interval I have constructed.” This statements is a direct probabilistic statement about the interval you have just created. Recall our key rule: it is forbidden in classical statistics to make direct probability statements about unobservable parameters. Memorize this. Your confidence interval only has an interpretation as part of an infinite set of other confidence intervals.
We have just hit upon the dirtiest open secret of classical statistics. There is no interpretation of your confidence interval other than this: the best you can say is that your interval either contains the actual value of μ or it does not, a statement which is a tautology, and, again therefore, true, but of no help (incidentally, Sullivan (2007) finally acknowledges this on p. 500). So what do you do with the interval you have just created? Why even bother, since it has no direct relation to the problem at hand? It’s even worse. Pick any two different numbers, say, 12 and 42. It is a true statement to say that this interval either contains μ or it does not for any statistical problem done by anybody with any data any time whatsoever (make sure you understand that before reading further).
The guy that invented confidence intervals, Dzerzij (Jerzy) Neyman, a statistician, knew about the interpretational problems of confidence intervals, and was concerned. But he was even more concerned about something called inductive arguments. An example due to Stove (1986): All the flames I have observed before have been hot (the premise); therefore, this flame will be hot (the conclusion). Neyman, and many other influential 20th century statisticians, rejected inductive arguments a basis for probability. They felt arguments like these were “groundless” or that inductive arguments were fallible because of the true statement that, for the flames example, there was nothing in the universe guaranteeing that this flame will be hot2. Inductive arguments are needed to make direct probabilistic statements about things like confidence intervals. If you reject them, then you cannot use probability. So Neyman, and those who followed him (which was nearly everybody), tried to take refuge in arguments like this: “Well, you cannot say that there is a 95% chance that the actual value of the parameter is in your interval; but if statisticians everywhere were to use confidence intervals, then in the long run, 95% of their intervals will contain their actual values.” Thirty-two extra credit points to those who can show the obvious flaw in this argument (see the homework).
The flaw in that argument was so obvious that it was evident to Neyman himself. And so, with nowhere else to turn, Neyman recommended a dodge and said this: “The statistician…may be recommended…to state that the value of the parameter μ is within (the just calculated interval)” merely by an act of will (Neyman (1937), quoted in Franklin (2001a)).
What you would like to be able to say is that “I have 95% (or whatever) confidence that this interval covers the true value of μ.” But you can never do this in classical statistics.
In R, to get the confidence interval of a normal distribution classically is a little more work than just getting the estimates, but it isn’t really that hard. This is for the appendicitis data, the White.Blood.Count (don’t forget to read the data in and attach it):
confint(glm(White.Blood.Count∼1))
The function confint calculates 95% confidence intervals. The inside function glm, with that funny argument ∼1, basically says, “The uncertainty in the variable should be quantified by a normal distribution.” Just take my word for it now; we’ll see this function later and this notation will become clear then. Anyway, after you run the command you will see something like this:
2.5 % 97.5 %
9.991874 10.818126
Ignore the word (Intercept), it is actually White.Blood.Count (this is because this function works for any variable name you care to enter). The 2.5 % and 97.5 % are like the quantiles; subtract 2.5 from 97.5 and get the length of the interval, which is 97.5%-2.5% = 95%.
We could use another R function and compute the confidence interval for σ, but it is not of great interest because later, we’ll see how to do all these things more or less automatically. Besides, we want to concentrate on what these intervals mean. If you’ve already forgotten, then go back and read this section from the beginning. One thing that is certain is that confidence intervals say nothing about the observables, the data x. If they say anything, they say something about the unobservable parameters. But what? The interval we computed for white blood count was about [10, 11]. This is an interval about estimated central parameter μ and not about the mean. We know the mean (it is…? find it in R). The confidence interval is an attempt to put a measure of precision on the guess μ. It says nothing about the mean, and nothing about actual values of white blood count. Never forget this.
5. Bayesian way
The idea behind modern statistics is that you quantify any and all uncertainty you have in anything using probability. We’ve already seen how to quantify uncertainty using probability for observables; that is, for actual data. That turns out to be done the same way classically and Bayesianly. This is what we did the first few Chapters, was it not? We wrote down some probability distribution, with known parameters, and made probability statements about observable data. Classical and Bayesian statistics begin to diverge when we start to talk about unknown parameters and how to make guesses about these parameters.
We made guesses classically by specifying some ad hoc function of the data, giving us θ; afterwards, we created a confidence interval for this guess. I stressed, heavily, that this confidence interval is not designed to express any actual uncertainty in θ, because that goes against the classical philosophy: which is that you cannot directly express uncertainty in unobservable parameters using probability.
In Bayesian statistics, you can, and must, express uncertainty in unobservable parameters using probability. How this works might sound complicated, and some of it is, but once you get how it works for, say, normal distributions, you will then know how it works for every other statistics problem in the world. This is not so for classical statistics, where you have to memorize a new set of ad hoc functions for every problem. In this way, Bayesian statistics is a vast simplification; however, before you can reach this simplification plateau, you initially have to climb up a steeper hill than you do classically. However, the good news is that there is only one hill to climb.
Let’s recall the normal probability distribution (density function):
(see the book)
written here as a function of x, or p(x|μ, σ, EN ) (we could have use N() as before; the actual letter does not matter). Do you remember probability rule number 4, or Bayes’s rule? If not, go back and re-read Chapter 2. Pay special attention to equation (6). I’ll wait here until you’re done.
Back? OK, let’s write equation (6) using different letters, so that
(see the book)
becomes
(see the book)
where B is now (μ, σ) and A is x. Remember, (μ, σ) is shorthand for the statement “The value of the central parameter is μ and the value of the variance parameter is σ”, and x is shorthand for the statement X = “The value of the observed data is x.” We already know how to write p(x|μ, σ, EN ) mathematically. Our goal is to discover how to write the left-hand side, which is the probability distribution of (μ, σ) given the data and EN . This quantifies our uncertainty in (μ, σ) given what we learned in the data (and considering the evidence EN ). In order to calculate the left-hand side, we then also need to know p(μ, σ|EN ). We also need p(x|EN ), but once we know p(μ, σ|EN ), it automatically pops out because of some math that need not concern us here.
What is p(μ, σ|EN )? Well, it quantifies our uncertainty in (μ, σ) before seeing any data, that is, it is only conditional on EN . p(μ, σ|EN ) is a probability distribution that you have to specify before you can get to the probability p(μ, σ|x, EN ). It also has an official name, which is the prior, because it’s what you know about (μ, σ) prior to adding in information in the data. Not surprisingly, then, p(μ, σ|x, EN ) is called the posterior, which is the probability distribution expressing everything we know, all our uncertainty, about (μ, σ) after having seen some data x.
How about the value of p(μ, σ|EN )? Well, it’s turns out to be a complicated situation, but the gist of it is that p(μ, σ|EN ) explains the probability of each possible value of (μ, σ), and since we initially know very little of (μ, σ), every possible value of (μ, σ) is more or less equally probable. This situation is called assigning a flat prior, the “flat” describing the shape of the probability distribution picture (i.e., a flat line)3. Once you have the prior, (There is more than one prior that you can use besides this “flat” one, but the differences it makes in the posteriors is minimal). Another problem is that the parameters are usually assumed to be continuous numbers, and if p(x|μ, σ, EN ), you can then calculate the posterior using equation (17). Technically, since we are saying (μ, σ) has a certain probability distribution, this is also information that we should keep note of, but we’ll append this on EN so that it now means “The uncertainty in the observable is quantified by a normal distribution and the prior on the parameters is ‘flat’.” If we need to be careful about this, and sometimes we do (not in this book), we can expand the notation to indicate the exact kind of prior we use.
Now here is another little secret: for very simple situations, the Bayesian results are the same as the classical results! No new calculations have to be learned or done!
After we take some old data, we can calculate our full uncertainty in (μ, σ) by drawing pictures of the probability distributions (we’ll do this later). If we are forced to pick just one “best” value, we would pick the arithmetic mean and standard deviation, exactly like in classical statistics. If we wanted to express our uncertainty a little more fully than just using one number (for each parameter), we could give the best number and an interval, some plus/minus bound on how certain that best value actually matches the true value of (μ, σ). Here is the best part: the confidence interval, which was meaningless before, is this interval, and is now called a credible interval. It has the natural interpretation that there is a 95% chance that the true value of the parameter lies in this interval. Isn’t that wild?
Before you start thinking, “Hey, if the results are the same, why did you go on and on and on about how confidence intervals are meaningless? All you did was to give them a new name! Big deal. You are wasting my time and trying to confuse me.” Hold on a minute, though. The Bayesian results are the same as the classical ones, but only for simple situations. The good news for you is, that in Stats 101, you hardly move beyond these very simple situations. Once you do move into the great statistical beyond, like using Binomial instead of normal distributions, the Bayesian methods really come into their own, and then you cannot you recall the discussion from Chapter 4, you know these can be a problem. We will ignore all these difficulties in this book. assume the classical computations give you the correct answer. I’ll talk about these techniques as we move along.
1I mean those people who were not formally trained in the mathematical subjects of probability and statistics. The vast numbers of people who compute statistics have not had this training beyond, say, a class given in a Psychology department by a professor who himself was not so trained, etc.
2To which you can argue; Ok, if you doubt it, stick your hand into this flame.
May 27th, 2008
It was one of those days yesterday. I got two chapters up, but did not give anybody a way to get them! Here it is
These are the last two “basics” Chapters. 6 first, and it is a little thin, so I’ll probably expand it later. It’s sort of a transition between probability where we know everything to statistics where we don’t. And by “everything” I mean the parameters of probability models. I want the reader to build up a little intuition before it starts to get rough.
The most important part of 6 is the homework, which I usually spend a lot of time with in class.
In a couple of days we start the good stuff. Book link.
CHAPTER 6
Normalities & Oddities
1. Standard Normal
Suppose x|m, s, EN ∼ N(m, s), then there turns out to be a trick that can make x easier to work with, especially if you have to do any calculations by hand (which, nowadays, will be rarely). Let z = (x-m)/s, then z|m, s, EN ∼ N(0, 1). It works for any m and s. Isn’t that nifty? Lots of fun facts about z can be found in any statistics textbook that weighs over 1 pound (these tidbits are usually in the form of impenetrable tables located in the back of the books).
What makes this useful is that Pr(z > 2|0, 1, EN ) ≈ Pr(z > 1.96|0, 1, EN ) = 0.025 and Pr(z < −2|0, 1, EN ) ≈ Pr(z < −1.96|0, 1, EN ) = 0.025: or, in words, the probability that z is bigger than 2 or less than negative 2 is about 0.05, which is a magic (I mean real voodoo) value in classical statistics. We already learned how to do this in R, last Chapter.
In Chapter 4, a homework question explained the rules of petanque, which is a game more people should play. Suppose the distance the boule lands from the cochonette is x centimeters. We do not know what x will be in advance, and so we (approximately) quantify our uncertainty in it using a normal distribution with parameters m = 0 cm and s = 10 cm. If x > 0 cm it means the boule lands beyond the cochonette, and if x < 0 cm is means the boule lands in front of the cochonette. You are out on the field playing, far from any computer, and the urge comes upon you to discover the probability that x > 30 cm. First thing to do is to calculate z which equals (30cm − 0cm)/10cm = 3 (the cm cancel). What is Pr(z > 3|0, 1, EN )? No idea; well, some idea. It must be less than 0.025, since we have all memorized that Pr(z > 2|0, 1, EN ) ≈ 0.025. The larger z is, the more improbable it becomes (right?). Let’s say as a guess 1%. When you get home, you can open R and plug in 1-pnorm(3) and see that the actually probability is 0.1%, so we were off by an order of magnitude (a power of 10), which is a lot, and which proves once again that computers are better at math than we are.
2. Nonstandard Normal
The standard normal example is useful for developing your probabilistic intuition. Since normal distributions are used so often, we will spend some more time thinking about some consequences of using them. Doing this will give you a better feel for how to quantify uncertainty.
Below is a picture of two normal distributions. The one with the solid line has m1 = 0 and s1 = 1; the dashed line has m2 = 0.5 and also s2 = 1. In other words, the two distributions differ only in their central parameter, they have the same variance parameter. Obviously, large values are more likely according to distribution 2, and smaller values are more likely given distribution 1, as a simple consequence of m2 > m1 . However, once we get to values of about x = 4 or so, it doesn’t look like the distributions are that different. (Cue the spooky music.) Or are they?.
Under the main picture are two others. The one on the left is exactly like the main picture, except that it focuses only on the range of x = 3.5 to x = 5. If we blow it up like this, we can see that it is still more likely to see large values of x using distribution 2.
How much more likely? The picture on the right divides the probabilities of seeing x or larger with distribution 2 by distribution 1, and so shows how much more likely it is to see larger values with distribution 2 than 1. For example, pick x = 4. It is about 7.5 times more likely to see an x = 4 or larger with distribution 2. That’s a lot! By the time we get out to x = 5, we are 12 times more likely to see values this large with distribution 2. The point is that even very small changes in the central parameters lead to large differences in the probabilities of “extreme”, values of x.
(see the book)
This next picture again shows two different distributions, this time with m1 = m2 = 0 with s1 = 1 and s1 = 1.1. In other words, both distributions have the same central parameters, but distribution 2 has a variance parameter that is slightly larger. The normal density plots do not look very different, do they? The dashed line, which is still distribution 2, has a peak slightly under distribution 1’s, but the differences looks pretty small.
(see the book)
The bottom panels are the same as before. The one on the left blows up the area where x > 3.5 and x < 5. A big difference still exists. And the ratio of probabilities is still very large. It's not shown, but the plot of the right would be duplicated (or mirrored, actually) if we looked at x > −5 and x < −3.5. It is more probable to see extreme events in either direction (positive or negative) using distribution 2.
The surprising consequence is that very small changes in either the central parameter or the variance parameter can lead to very large differences at the extremes. Examples of these phenomena are easily found in real life, but my heightened political sensitivity precludes me from publicly pointing any of these out.
3. Intuition
We have learned probability and some formal distributions, but we have not yet moved to statistics. Before we do so, let us try to develop some intuition about the kinds of problems and solutions we will see before getting to technicalities. There are a number of concepts that will be important, but I don’t want to give them a name, because there is no need to memorize jargon, while it is incredibly important that you develop a solid under- standing of uncertainty.
The well-known Uncle Ted Nugent’s chain of Kill ‘em and Grill ‘em Vension Burger restaurants sell both Coke and Pepsi, and their internal audit shows they sell about an equal amount of each. The busy Times Square branch of the chain has about 5000 customers a day, while the store in tiny Gaylord, Michigan sees only about 100 customers. Which location is more likely to sell, on any given day, at least 2 times more Pepsi than Coke?
A useful technique for solving questions like this is exaggeration. For instance, the question is asking about a difference in location. What differs between those places? Only one thing, the number of customers. One site gets about 5000 people a day, the other only 100. Let’s exaggerate that difference and solve a simpler problem. For example, suppose Times Square still gets 5000 a day, but Gaylord only gets 1 a day. The information is that selling a Coke is roughly equal to the probability of selling a Pepsi. This means that, at Gaylord, to that 1 customer on that day, they will either sell 1 Coke or 1 Pepsi. If they sell a Pepsi, Gaylord has certainly sold more than 2 times as much Pepsi as Coke. The chance of that happening is 50%. What is two times as much Pepsi as Coke at Times Square? A lot more Pepsi, certainly. So it’s far more likely for Gaylord to sell a greater proportion of Pepsi because they see fewer customers. The lesson is that when the “sample size” is small, we are more likely to see extreme events.
What is the length of the first Chinese Emperor Qin Shi Huangdi’s nose? You don’t know? Well, you can make a guess. How likely is it that your guess is correct? Not very likely. Suppose that you decide to ask everybody you know to also guess, and then average all the answers together in an attempt to get a better guess. How likely is it that this averaged-guess is perfectly correct? No more likely. If you haven’t a clue about the nose, and nobody else does either, than averaging ignorance is no better than single ignorance. The lesson is that just because a large group of people agree on an opinion, it is not necessarily more probable that that opinion, or average of opinions, is correct. Uninformed opinion of a large group of people is not necessarily more likely to be correct than the opinion of the lone nut job on the corner. Think about this the next time you hear the results of a poll or survey.
You already posses other probabilistic intuition. For example, suppose, given some evidence E, the probability of A is 0.0000001 (A is something that might be given many opportunities to happen, e.g. winning the lottery). How often will A happen? Right. Not very often. But if you give A a lot of chances to occur, will A eventually happen? It’s very likely to.
Every player in petanque gets to throw three boules. What are the chances that I get all three within 5 cm? This is a compound problem, so let’s break it apart. How do we find out how likely it is to be within 5 cm of the cochonette? Well, that means the boule can be 5 cm in front of the cochonette, right near it, or up to 5cm beyond it. The chance of this happening is Pr(−5cm < x < 5cm|m = 0cm, s = 10cm, EN ). We learned how to calculate the probability of being in an interval last chapter:
pnorm(5,0,10)-pnorm(-5,0,10).
This equals about 0.38, which is the chance that one boule lands within, or +/- 5 cm, from the cochonette. What is the chance that all of them land that close? Well, that means the first one does and the second one and the third. What probability rule do we use now? The second, which tells us to multiple the probabilities together, which is 0.383 ≈ 0.14. The important thing to recall, when confronted with problems of this sort: do not panic. Try to break apart the complex problem into bite-size pieces.
May 22nd, 2008
Update: 21 May 4:45 am. I forgot to actually upload the file until right this moment. Thanks to Mike and Harry for the reminder.
Chapter 5 is ready to go.
This is purely a mechanical chapter, introducing R. Thrilling reading, it is not. But it’s necessary to learn in order to be able to carry out the analysis in later chapters. The book website is not fully up; only the datasets are there. To learn to install R, just look on the R website.
I’ll be posting Chapters 6 and 7 in short order and then we finally get to the good stuff.
R
1. R
R is a fantastic, hugely supported, rapidly growing, infinitely extensible, operating-system agnostic, free and open source statistical software platform. Nearly everybody who is anybody uses R, and since I want you to be somebody, you will use it, too. Some things in R are incredibly easy to do; other tasks are bizarrely difficult. Most of what makes R hard for the beginner is the same stuff that makes any piece of software hard; that is, getting used to expressing your statistical desires in computerese. As such an environment can be strange and perplexing at first, some students experience a kind of peculiar stress that is best described by example. Here is a video from a Germany showing a young statistics student who experienced trouble understanding R:
http://youtube.com/watch?v=PbcctWbC8Q0
Be sure that this doesn’t happen to you. Remember what Douglas Adams said: Don’t panic.
The best way to start is by going to r-project.org and click the CRAN under the Download heading. You can’t miss it. After that, you have to choose a mirror, which means one of the hundreds of computers around the world that host the software. Obviously, pick a site near you. Once that’s done, and choose your platform (your operating system, like Linux or one of the others), and then choose the base package. Step-by-step instructions are at this book’s website: wmbriggs.com/book. It is no more difficult to install than any other piece of software.
This is not the place to go over all the possibilities of R; just the briefest introduction will be given, because there are far better places available online (see the book website for links). But there are a few essential commands that you should not do without.
These are
| Command |
Description |
| help(command) |
Does the obvious: always scroll down to
the bottom of the help to see examples of
the command.
|
| ?command |
Same as help() |
| apropos(’string’) |
If you cannot remember the name of a
command—and I always forget—but re-
member is started with co–something,
then just type apropos(’co’) and you’ll
get a complete list of commands that have
co anywhere in their names. |
| c() |
This is the concatenation function: typing
c(1,2) concatenates a 2 to 1, or sticks on
the end 1 the number 2, so that we have
a vector of numbers. |
The Appendix gives a fuller list of R commands.
It is important to understand that R is a command-line language, which we may interpret as meaning that all commands in R are functions which must be typed into the console. These are objects that are a command name plus a left and right parenthesis, with variables (called arguments) stuck in between, thus: plot(x,y). Remember that you are dealing with computers, which are literal, intolerant creatures (much like the people who want to ban smoking), and so cannot abide even the slightest deviation from its expectations. That means, if instead of plot(x,y), you type lot(x,y), or plot x,y), or plot(,y), or plot(x,y things will go awry. R will try to give you an idea of what went wrong by giving you an error message. Except in cases like that last typo, which will cause you to develop stress lines, because all you’ll see is this
+
and every attempt of yours to type anything new, or hit enter 100 times, will not do a thing except give you more lines of + or other screwy errors. Because why? Because you typed plot(x,y; that is, you typed a left parenthesis (right before the x) and you never “closed” it with a right parenthesis, and R will simply wait forever for you to type one in.
The solution is to enter a right parenthesis, or hit
ctrl+c
which means the control key plus the c key simultaneously, which “breaks” the current computation.
Using R means that you have to memorize (!) and type in commands instead of using a graphical user interface (GUI), which is the standard point-and-click screen with which you are probably familiar. It is my experience that students who are not used to computers start freaking out at this point; however, there is no need to. I have made everything very, very easy and all you have to do is copy what you see in the book to the R screen. All will be well. I promise.
GUIs are very nice things, incidentally, and R has one that you can download and play with. It is called the R Commander. Like all GUIs, some very basic functionality is included that allows you to, well, point and click and get a result. Problem is, the very second you want to do something different than what is available from the GUI, you are stuck. With statistics, we often want to do something differently, so we will stick with the command line.
2. R binomially
By now, you are eagerly asking yourself:”“Can R help up with those binomial calculations like in the Thanksgiving example?” Let’s type apropos(’bino’) and see, because, after all, ‘bino’ is something like binomial. The most likely function is called binomial, so let’s type ?binomial and see how it works. Uh oh. Weird words about “family objects” and the function glm(), and that doesn’t sound right. What about one of the functions like dbinom()? Jackpot. We’ll look at these in detail, since it turns out that this structure of four functions is the same for every distribution. The functions are in this table:
|
dbinom |
The probability of density function: given the size,
or n, and prop, or p, this calculates the probability
that we see x successes; this is equation (11). |
| pbinom |
The distribution function, which calculates the probability that the number of successes is less than or
equal to some a. |
| qbinom |
This is the “quantile” function, which calculates,
given a probability from the distribution function,
which value of q it is associated with. This will be
made clear with some examples with the normal distribution later. |
| rbinom |
This generates a “random” binomial number; and
since random means unknown, this means it gener-
ates a number that is unknown in some sense; we’ll
talk about this later. |
Let’s go back to the Thanksgiving example, which used a binomial. Moe can calculate, given n = size = 3, p = prob = 0.1,
his probabilities using R:
dbinom(0,3,.1)
which gives the probability of taking nobody along for the ride. The answer is [1] 0.729. The “[1] in front of the number just means that you are only looking at line number 1. If you asked for dozens of probabilities, for example, R would space them out over several lines. Let’s now calculate the probability of taking just 0, just 1, etc.
dbinom(c(0,1,2,3),3,.1)
where we have “nested” two functions into one: the first is the concatenation function c(), where we have stuck the numbers 0 through 3 together, and which shows you the dbinom() function can calculate more than one probability at a time. What pops out is
[1] 0.729 0.243 0.027 0.001;
that is, the exact values we got above for taking 0 or 1 or 2 etc. along for the ride. Now we can look at the distribution function:
pbinom(c(0,1,2,3),3,.1);
and we get
[1] 0.729 0.972 0.999 1.000.
This is the probability of taking 0 or less, 1 or less, 2 or less, and 3 or less. The last probability very obviously has to be 1, and will always be 1 for any binomial (as long as the last value in the function c(0,1,…,n) equals n).
There turns out to be a shortcut to typing the concatenation function for simple numbers, and here it is:
c(0,1,2,…,n) = 0:n.
So we can rewrite the first function as dbinom(0:3,3,.1) and get the same results.
We can nest functions again and make pretty pictures
plot(dbinom(0:3,3,.1))
And that’s it for any binomial function. Isn’t that simple? (The answer is yes.) The commands never change for any binomial you want to do.
3. R normally
Can R do normal distributions as well? Can it! Let’s type in apropos(’normal’) and see what we get. A lot of gibberish, that’s what. Where’s the normal distribution? Well, it turns out that computer programmers are a lazy bunch, and they often do not use all the letters of a word to name a function (too much typing). Let’s try apropos(’norm’) instead (which no matter what should give us at least as many results, right? This is a question of logic, not computers.). Bingo. Among all the rest, we see dnorm and pnorm etc., just like with the biomial. Now type ?dnorm so we can learn about our fun function. Same layout as the binomial; only difference being we need to supply a “mean” and “sd” (the m and s). Sigh. This is an example of R being naughty and misusing the terminology that I earlier forbade: m and s are not a mean and standard deviation. It’s a trap too many fall into. We’ll work with it, but just remember “mean” and “sd” actually imply our parameters m and s.
You will recall from our discussion of normals that we cannot compute a probability of seeing a single number (and if you don’t remember, shame on you: go back and read Chapter 4). The function dnorm does not give you this number, because that probability is always 0; instead, it gives you a “density”, which means little to us. But we can calculate the probability of values being in some interval using the pnorm function. For example, to calculate Pr(x < 10|m = 10, s = 20, EN ), use
pnorm(10,10,20)
and you should see [1] 0.5. But you already knew that would be the answer before you typed it in, right? (Right?) Let’s try a trickier one: Pr(x < 0|m = 10, s = 20, EN ); type pnorm(0,10,20) and get [1] 0.3085375. So what is this probability: Pr(x > 0|m = 10, s = 20, EN ) (x greater than 0)? Think about it. x can either be less than or greater than 0; the probability it is so is 1. So Pr(x < 0|m = 10, s = 20, EN ) + Pr(x > 0|m = 10, s = 20, EN ) = 1. Thus, Pr(x < 0|m = 10, s = 20, EN ) Pr(x > 0|m = 10, s = 20, EN ) = 1 − Pr(x < 0|m = 10, s = 20, EN ). We can get that in R by typing
1-pnorm(0,10,20)
and you should get [1] 0.6914625, which is 1 − 0.3085375 as expected.
By the way, if you are starting to feel the onset of a freak out, and wonder "Why, O why, can't he give us a point and click way to do this!" Because, dear reader, a point and click way to do this does not exist. Stop worrying so much. You'll get it.
What is Pr(15 < x < 18|m = 15, s = 5, EN ) (which might be reasonable numbers for the temperature example)? Any interval splits the data into three parts: the part less than the lower bound (15), the part of the interval itself (15-18), and the part larger than the upper bound (18). We already know how to get Pr(x < 15|m = 15, s = 5, EN ), which is pnorm(15,15,5), and which equals 0.5. We also know how to get Pr(x > 18|m = 15, s = 5, EN ), which is 1-pnorm(18,15,5), and which equals 0.2742531. This means that Pr(x < 15 or x > 18|m = 15, s = 5, EN ), using probability rule number 1, is 0.5 + 0.2742531 = 0.7742531. Finally, 0.7742531 is the probability of not being in the interval, so the probability of being in the interval must be one minus this, or 1 − 0.7742531 = 0.2257469. A lot of work. We could have jumped right to it by typing
pnorm(18,15,5)-pnorm(15,15,5).
This is the way you write the code to compute the probability of any interval—remembering to input your own m and s of course!
4. Advanced
. You don’t need to do this section, because it is somewhat more complicated. Not much, really, but enough that you have to think more about the computer than you do the probability.
Our goal is to plot the picture of a normal density. The function dnorm(x,15,5) will give you the value of the normal density, with an m = 15 and s = 5, for some value of x. To picture the normal, which is a picture of densities for a range of x, we somehow have to specify this range. Unfortunately, there is no way to know in advance which range you want to plot, so getting the exact picture you want takes some work. Here is one way:
x = seq(-4,4,.01)
which gives us a sequence of numbers from -4 to 4 in increments of 0.01. Thus, x = −4.00, −3.99, −3.99, . . . , 4. Calculating the density of each of these values of x is easy:
dnorm(x)
where you will have noticed that I did not type a m or s. Type ?dnorm again. It reads dnorm(x, mean=0, sd=1, log = FALSE). Ignoring the log = FALSE bit, we can see that R supplies helpfully default values of the parameters. They are default, because if you are happy with the values chosen, you do not have to type in your own. In this case, m = 0 and s = 1, which is called a standard normal. Anyway, to get the plot is now easy:
plot(x,dnorm(x),type=’l’)
This means, for every value of x, plot the value of dnorm at that value. I also changed the plot type to a line with type=’l’ , and which makes the graph prettier. Try doing the plot without this argument and see what you get.
May 20th, 2008
Chapter 4 is ready to go.
This is where it starts to get weird. The first part of the chapter introduces the standard notation of “random” variables, and then works through a binomial example, which is simple enough.
Then come the so-called normals. However, they are anything but. For probably most people, it will be the first time that they hear about the strange creatures called continuous numbers. It will be more surprising to learn that not all mathematicians like these things or agree with their necessity, particularly in problems like quantifying probability for real observable things.
I use the word “real” in its everyday, English sense of something that is tangible or that exists. This is because mathematicians have co-opted the word “real” to mean “continuous”, which in an infinite amount of cases means “not real” or “not tangible” or even “not