Probably Overthinking It

Self-organized criticality and holistic models

2012-02-24T08:41:00.001-08:00

Long-tailed distributions of some physical quantities: for example, in freezing water the distribution of crystal sizes is characterized by a power law.
Fractal geometries: freezing water tends to form fractal patterns---the canonical example is a snowflake. Fractals are characterized by self-similarity; that is, parts of the pattern resemble scaled copies of the whole.
Variations in time that exhibit pink noise: what we call ``noise'' is a time series with many frequency components. In ``white'' noise, all of the components have equal power. In ``pink'' noise, low-frequency components have more power than high-frequency components. Specifically, the power at frequency f is proportional to 1/f. Visible light with this power spectrum looks pink, hence the name.

Critical systems are usually unstable. For example, to keep water in a partially frozen state requires active control of the temperature. If the system is near the critical temperature, a small deviation tends to move the system into one phase or the other.

Many natural systems exhibit characteristic behaviors of criticality, but if critical points are unstable, they should not be common in nature. This is the puzzle Bak, Tang and Wiesenfeld address. Their solution is called self-organized criticality (SOC), where ``self-organized'' means that from any initial condition, the system tends to move toward a critical state, and stay there, without external control.

As an example, they propose a model of a sand pile. The model is not realistic, but it has become the standard example of self-organized criticality.

The model is a 2-D cellular automaton where the state of each cell represents the slope of a part of a sand pile. During each time step, each cell is checked to see whether it exceeds some critical value. If so, an ``avalanche'' occurs that transfers sand to neighboring cells; specifically, the cell's slope is decreased by 4, and each of the 4 neighbors is increased by 1. At the perimeter of the grid, all cells are kept at zero slope, so (in some sense) the excess spills over the edge.

Bak et al. let the system run until it is stable, then observe the effect of small perturbations; they choose a cell at random, increment its value by 1, and evolve the system, again, until it stabilizes.

For each perturbation, they measure the total number of cells that are affected by the resulting avalanche. Most of the time it is small, usually 1. But occasionally a large avalanche affects a substantial fraction
of the grid. The distribution of turns out to be long-tailed, which supports the claim that the system is in a critical state.

[Think Complexity presents the details of this model and tests for long-tailed distributions, fractal geometry, and 1/f noise. For this excerpt, I'll skip to the discussion at the end of the chapter.]

Reductionism and Holism

The original paper by Bak, Tang and Wiesenfeld is one of the most frequently-cited papers in the last few decades. Many new systems have been shown to be self-organized critical, and the sand-pile model, in particular, has been studied in detail.

As it turns out, the sand-pile model is not a very good model of a sand pile. Sand is dense and not very sticky, so momentum has a non-negligible effect on the behavior of avalanches. As a result, there are fewer very large and very small avalanches than the model predicts, and the distribution is not long tailed.

Bak has suggested that this observation misses the point. The sand pile model is not meant to be a realistic model of a sand pile; it is meant to be a simple example of a broad category of models.

To understand this point, it is useful to think about two kinds of models, reductionist and holistic. A reductionist model describes a system by describing its parts and their interactions. When a reductionist model is used as an explanation, it depends on an analogy between the components of the model and the components of the system.

For example, to explain why the ideal gas law holds, we can model the molecules that make up a gas with point masses, and model their interactions as elastic collisions. If you simulate or analyze this model, you find that it obeys the ideal gas law. This model is satisfactory to the degree that molecules in a gas behave like molecules in the model. The analogy is between the parts of the system and the parts of the model.

Holistic models are more focused on similarities between systems and less interested in analogous parts. A holistic approach to modeling often consists of two steps, not necessarily in this order:

1. Identify a kind of behavior that appears in a variety of systems.

2. Find the simplest model that demonstrates that behavior.

For example, in The Selfish Gene, Richard Dawkins suggests that genetic evolution is just one example of an evolutionary system. He identifies the essential elements of the category---discrete replicators, variability and differential reproduction---and proposes that any system that has these elements displays similar behavior, including complexity without design. As another example of an evolutionary system, he proposes memes, which are thoughts or behaviors that are ``replicated'' by transmission from person to person. As memes compete for the resource of human attention, they evolve in ways that are similar to genetic evolution.

Critics of memetics have pointed out that memes are a poor analogy for genes. Memes differ from genes in many obvious ways. But Dawkins has argued that these differences are beside the point because memes are not supposed to be analogous to genes. Rather, memetics and genetics are examples of the same category---evolutionary systems. The differences between them emphasize the real point, which is that evolution is a general model that applies to many seemingly disparate systems. The logical structure of this argument is shown in this diagram:

Bak has made a similar argument that self-organized criticality is a general model for a broad category of systems. According to Wikipedia, ``SOC is typically observed in slowly-driven non-equilibrium systems with extended degrees of freedom and a high level of nonlinearity.''

Many natural systems demonstrate behaviors characteristic of critical systems. Bak's explanation for this prevalence is that these systems are examples of the broad category of self-organized criticality. There are two ways to support this argument. One is to build a realistic model of a particular system and show that the model exhibits SOC. The second is to show that SOC is a feature of many diverse models, and to identify the essential characteristics those models have in common.

The first approach, which I characterize as reductionist, can explain the behavior of a particular system. The second, holistic, approach, explains the prevalence of criticality in natural systems. They are different models with different purposes.

For reductionist models, realism is the primary virtue, and simplicity is secondary. For holistic models, it is the other way around.

I am using "reductionism" and "holism" here is a descriptive sense, not as technical labels for these models. For more general discussion of these terms, see http://en.wikipedia.org/wiki/Reductionism and http://en.wikipedia.org/wiki/Holism.

SOC, causation and prediction

If a stock market index drops by a fraction of a percent in a day, there is no need for an explanation. But if it drops 10, people want to know why. Pundits on television are willing to offer explanations, but the real answer may be that there is no explanation.

Day-to-day variability in the stock market shows evidence of criticality: the distribution of value changes is long-tailed and the time series exhibits noise. If the stock market is a self-organized critical system, we should expect occasional large changes as part of the ordinary behavior of the market.

The distribution of earthquake sizes is also long-tailed, and there are simple models of the dynamics of geological faults that might explain this behavior. If these models are right, they imply that large earthquakes are unexceptional; that is, they do not require explanation any more than small earthquakes do.

Similarly, Charles Perrow has suggested that failures in large engineered systems, like nuclear power plants, are like avalanches in the sand pile model. Most failures are small, isolated and harmless, but occasionally a coincidence of bad fortune yields a catastrophe. When big accidents occur, investigators go looking for the cause, but if Perrow's ``normal accident theory'' is correct, there may be no cause.

These conclusions are not comforting. Among other things, they imply that large earthquakes and some kinds of accidents are fundamentally unpredictable. It is impossible to look at the state of a critical system and say whether a large avalanche is ``due.'' If the system is in a critical state, then a large avalanche is always possible. It just depends on the next grain of sand.

In a sand-pile model, what is the cause of a large avalanche? Philosophers sometimes distinguish the proximate cause, which is most immediately responsible, from the ultimate cause, which is, for whatever reason, considered the true cause.

In the sand-pile model, the proximate cause of an avalanche is a grain of sand, but the grain that causes a large avalanche is identical to any other grain, so it offers no special explanation. The ultimate cause of a large avalanche is the structure and dynamics of the systems as a whole: large avalanches occur because they are a property of the system.

Many social phenomena, including wars, revolutions, epidemics, inventions and terrorist attacks, are characterized by long-tailed distributions. If the reason for these distributions is that social systems are critical, that suggests that major historical events may be fundamentally unpredictable and unexplainable.

Questions

[Think Complexity can be used as a textbook, so it includes exercises and topics for class discussion. Here are some ideas for discussion and further reading.]

1. In a 1996 paper in Nature, Frette et al report the results of experiments with rice piles (http://www.nature.com/nature/journal/v379/n6560/abs/379049a0.html). They find that some kinds of rice yield evidence of critical behavior, but others do not.

Similarly, Pruessner and Jensen studied large-scale versions of the forest fire model (using an algorithm similar to Newman and Ziff's). In their 2004 paper, ``Efficient algorithm for the forest fire model,'' they present evidence that the system is not critical after all (http://pre.aps.org/abstract/PRE/v70/i6/e066707). How do these results bear on Bak's claim that SOC explains the prevalence of critical phenomena in nature?

2. In The Fractal Geometry of Nature, Benoit Mandelbrot proposes what he calls a ``heretical'' explanation for the prevalence of long-tailed distributions in natural systems (page 344). It may not be, as Bak suggests, that many systems can generate this behavior in isolation. Instead there may be only a few, but there may be interactions between systems that cause the behavior to propagate.

To support this argument, Mandelbrot points out:

The distribution of observed data is often ``the joint effect of a fixed underlying 'true distribution' and a highly variable 'filter.'''
Long-tailed distributions are robust to filtering; that is, ``a wide variety of filters leave their asymptotic behavior unchanged.''

What do you think of this argument? Would you characterize it as reductionist or holist?

3. Read about the ``Great Man'' theory of history at http://en.wikipedia.org/wiki/Great_man_theory. What implication does self-organized criticality have for this theory?

Think Complexity: Part Four

2012-02-14T08:55:00.000-08:00

My new book, Think Complexity, will be published by O'Reilly Media in March. For people who can't stand to wait that long, I am publishing excerpts here. If you really can't wait, you can read the free version at thinkcomplex.com.

And we need a blurb. Think Complexity goes to press soon and we have a space on the back cover for a couple of endorsements. If you like the book and have something quotable to say about it, let me know. Thanks!

In Part One I outline the topics in Think Complexity and contrasted a classical physical model of planetary orbits with an example from complexity science: Schelling's model of racial segregation.

In Part Two I outline some of the ways complexity differs from classical science. In Part Three, I describe differences in the ways complex models are used, and their effects in engineering and (of all things) epistemology.

In this installment, I pull together discussions from two chapters: the Watts-Strogatz model of small world graphs, and the Barabasi-Albert model of scale free networks. But it all starts with Stanley Milgram.

Stanley Milgram

Stanley Milgram was an American social psychologist who conducted two of the most famous experiments in social science, the Milgram experiment, which studied people's obedience to authority (http://en.wikipedia.org/wiki/Milgram_experiment) and the Small World Experiment (http://en.wikipedia.org/wiki/Small_world_phenomenon), which studied the structure of social networks.

In the Small World Experiment, Milgram sent a package to several randomly-chosen people in Wichita, Kansas, with instructions asking them to forward an enclosed letter to a target person, identified by name and occupation, in Sharon, Massachusetts (which is the town near Boston where I grew up). The subjects were told that they could mail the letter directly to the target person only if they knew him personally; otherwise they were instructed to send it, and the same instructions, to a relative or friend they thought would be more likely to know the target person.

Many of the letters were never delivered, but of the ones that were it turned out that the average path length---the number of times the letters were forwarded---was about six. This result was taken to confirm previous observations (and speculations) that the typical distance between any two people in a social network is about ``six degrees of separation.''

This conclusion is surprising because most people expect social networks to be localized---people tend to live near their friends---and in a graph with local connections, path lengths tend to increase in proportion to geographical distance. For example, most of my friends live nearby, so I would guess that the average distance between nodes in a social network is about 50 miles. Wichita is about 1600 miles from Boston, so if Milgram's letters traversed typical links in the social network, they should have taken 32 hops, not six.

Watts and Strogatz

In 1998 Duncan Watts and Steven Strogatz published a paper in Nature, ``Collective dynamics of 'small-world' networks,'' that proposed an explanation for the small world phenomenon. You can download it from http://www.nature.com/nature/journal/v393/n6684/abs/393440a0.html.

Watts and Strogatz started with two kinds of graph that were well understood: random graphs and regular graphs. They looked at two properties of these graphs, clustering and path length.

Clustering is a measure of the ``cliquishness'' of the graph. In a graph, a clique is a subset of nodes that are all connected to each other; in a social network, a clique is a set of friends who all know each other. Watts and Strogatz defined a clustering coefficient that quantifies the likelihood that two nodes that are connected to the same node are also connected to each other.

Path length is a measure of the average distance between two nodes, which corresponds to the degrees of separation in a social network.

Their initial result is what you might expect: regular graphs have high clustering and high path lengths; random graphs with the same size tend to have low clustering and low path lengths. So neither of these is a good model of social networks, which seem to combine high clustering with short path lengths.

Their goal was to create a generative model of a social network. A generative model tries to explain a phenomenon by modeling the process that builds or leads to the phenomenon. In this case Watts and Strogatz proposed a process for building small-world graphs:

Start with a regular graph with n nodes and degree k. Watts and Strogatz start with a ring lattice, which is a kind of regular graph. You could replicate their experiment or try instead a graph that is regular but not a ring lattice.
Choose a subset of the edges in the graph and ``rewire'' them by replacing them with random edges. Again, you could replicate the procedure described in the paper or experiment with alternatives. The proportion of edges that are rewired is a parameter, p, that controls how random the graph is. With p=0, the graph is regular; with p=1 it is random.

Watts and Strogatz found that small values of p yield graphs with high clustering, like a regular graph, and low path lengths, like a random graph.

Barabasi and Albert

In 1999 Barabasi and Albert published a paper in Science, ``Emergence of Scaling in Random Networks,'' that characterizes the structure (also called ``topology'') of several real-world networks, including graphs that represent the interconnectivity of movie actors, world-wide web (WWW) pages, and elements in the electrical power grid in the western United States. You can download the paper from http://www.sciencemag.org/content/286/5439/509.

They measure the degree (number of connections) of each node and compute P(k), the probability that a vertex has degree k; then they plot P(k) versus k on a log-log scale. The tail of the plot fits a straight line, so they conclude that it obeys a power law; that is, as k gets large, P(k) is asymptotic to k^(-γ), where γ is a parameter that determines the rate of decay.

They also propose a model that generates random graphs with the same property. The essential features of the model, which distinguish it from the model and the Watts-Strogatz model, are:

Growth: Instead of starting with a fixed number of vertices, Barabasi and Albert start with a small graph and add vertices gradually.

Preferential attachment: When a new edge is created, it is more likely to connect to a vertex that already has a large number of edges. This ``rich get richer'' effect is characteristic of the growth patterns of some real-world networks.

Finally, they show that graphs generated by this model have a distribution of degrees that obeys a power law. Graphs that have this property are sometimes called scale-free networks; see http://en.wikipedia.org/wiki/Scale-free_network. That name can be confusing because it is the distribution of degrees that is scale-free, not the network.

In order to maximize confusion, distributions that obey a power law are sometimes called scaling distributions because they are invariant under a change of scale. That means that if you change the units the quantities are expressed in, the slope parameter, γ, doesn't change. You can read http://en.wikipedia.org/wiki/Power_law for the details, but it is not important for what we are doing here.

Explanatory models

We started the discussion of networks with Milgram's Small World Experiment, which shows that path lengths in social networks are surprisingly small; hence, ``six degrees of separation''. When we see something surprising, it is natural to ask ``Why?'' but sometimes it's not clear what kind of answer we are looking for.

One kind of answer is an explanatory model. The logical structure of an explanatory model is:

In a system, S, we see something observable, O, that warrants explanation.
We construct a model, M, that is analogous to the system; that is, there is a correspondence between the elements of the model and the elements of the system.
By simulation or mathematical derivation, we show that the model exhibits a behavior, B, that is analogous to O.
We conclude that S exhibits O because S is similar to M, M exhibits B, and B is similar to O.

At its core, this is an argument by analogy, which says that if two things are similar in some ways, they are likely to be similar in other ways. Argument by analogy can be useful, and explanatory models can be satisfying, but they do not constitute a proof in the mathematical sense of the word.

Remember that all models leave out, or ``abstract away'' details that we think are unimportant. For any system there are many possible models that include or ignore different features. And there might be models that exhibit different behaviors, B, B' and B'', that are similar to O in different ways. In that case, which model explains O?

The small world phenomenon is an example: the Watts-Strogatz (WS) model and the (BA) model both exhibit small world behavior, but they offer different explanations:

The WS model suggests that social networks are ``small'' because they include both strongly-connected clusters and ``weak ties'' that connect clusters.
The BA model suggests that social networks are small because they include nodes with high degree that act as hubs, and that hubs grow, over time, due to preferential attachment.

As is often the case in young areas of science, the problem is not that we have no explanations, but too many.

Questions

[Think Complexity can be used as a textbook, so it includes exercises and topics for class discussion. Here are some ideas for discussion and further reading.]

Are these explanations compatible; that is, can they both be right? Which do you find more satisfying as an explanation, and why? Is there data you could collect, or an experiment you could perform, that would provide evidence in favor of one model over the other?

Choosing among competing models is the topic of Thomas Kuhn's essay, ``Objectivity, Value Judgment, and Theory Choice.'' You can download it here in PDF. What criteria does Kuhn propose for choosing among competing models? Do these criteria influence your opinion about the WS and BA models? Are there other criteria you think should be considered?

Think Complexity: Part Three

2012-02-06T08:23:00.000-08:00

My new book, Think Complexity, will be published by O'Reilly Media in March. For people who can't stand to wait that long, I am publishing excerpts here. If you really can't wait, you can read the free version at thinkcomplex.com.

In Part One I outline the topics in Think Complexity and contrasted a classical physical model of planetary orbits with an example from complexity science: Schelling's model of racial segregation.

In Part Two I outline some of the ways complexity differs from classical science. In this installment, I describe differences in the ways complex models are used, and their effects in engineering and (of all things) epistemology.

A new kind of model

Complex models are often appropriate for different purposes and interpretations:

Predictive→explanatory: Schelling's model of segregation might shed light on a complex social phenomenon, but it is not useful for prediction. On the other hand, a simple model of celestial mechanics can predict solar eclipses, down to the second, years in the future.

Realism→instrumentalism: Classical models lend themselves to a realist interpretation; for example, most people accept that electrons are real things that exist. Instrumentalism is the view that models can be useful even if the entities they postulate don't exist. George Box wrote what might be the motto of instrumentalism: ``All models are wrong, but some are useful."

Reductionism→holism: Reductionism is the view that the behavior of a system can be explained by understanding its components. For example, the periodic table of the elements is a triumph of reductionism, because it explains the chemical behavior of elements with a simple model of the electrons in an atom. Holism is the view that some phenomena that appear at the system level do not exist at the level of components, and cannot be explained in component-level terms.

A new kind of engineering

I have been talking about complex systems in the context of science, but complexity is also a cause, and effect, of changes in engineering and the organization of social systems:

Centralized→decentralized: Centralized systems are conceptually simple and easier to analyze, but decentralized systems can be more robust. For example, in the World Wide Web clients send requests to centralized servers; if the servers are down, the service is unavailable. In peer-to-peer networks, every node is both a client and a server. To take down the service, you have to take down every node.

Isolation→interaction: In classical engineering, the complexity of large systems is managed by isolating components and minimizing interactions. This is still an important engineering principle; nevertheless, the availability of cheap computation makes it increasingly feasible to design systems with complex interactions between components.

One-to-many→many-to-many: In many communication systems, broadcast services are being augmented, and sometimes replaced, by services that allow users to communicate with each other and create, share, and modify content.

Top-down→bottom-up: In social, political and economic systems, many activities that would normally be centrally organized now operate as grassroots movements. Even armies, which are the canonical example of hierarchical structure, are moving toward devolved command and control.

Analysis→computation: In classical engineering, the space of feasible designs is limited by our capability for analysis. For example, designing the Eiffel Tower was possible because Gustave Eiffel developed novel analytic techniques, in particular for dealing with wind load. Now tools for computer-aided design and analysis make it possible to build almost anything that can be imagined. Frank Gehry's Guggenheim Museum Bilbao is my favorite example.

Design→search: Engineering is sometimes described as a search for solutions in a landscape of possible designs. Increasingly, the search process can be automated. For example, genetic algorithms explore large design spaces and discover solutions human engineers would not imagine (or like). The ultimate genetic algorithm, evolution, notoriously generates designs that violate the rules of human engineering.

A new kind of thinking

We are getting farther afield now, but the shifts I am postulating in the criteria of scientific modeling are related to 20th Century developments in logic and epistemology.

Aristotelian logic→many-valued logic: In traditional logic, any proposition is either true or false. This system lends itself to math-like proofs, but fails (in dramatic ways) for many real-world applications. Alternatives include many-valued logic, fuzzy logic, and other systems designed to handle indeterminacy, vagueness, and uncertainty. Bart Kosko discusses some of these systems in Fuzzy Thinking.

Frequentist probability →Bayesianism: Bayesian probability has been around for centuries, but was not widely used until recently, facilitated by the availability of cheap computation and the reluctant acceptance of subjectivity in probabilistic claims. Sharon Bertsch McGrayne presents this history in The Theory That Would Not Die.

Objective→ subjective: The Enlightenment, and philosophic modernism, are based on belief in objective truth; that is, truths that are independent of the people that hold them. 20th Century developments including quantum mechanics, Godel's Incompleteness Theorem, and Kuhn's study of the history of science called attention to seemingly unavoidable subjectivity in even ``hard sciences'' and mathematics. Rebecca Goldstein presents the historical context of Godel's proof in Incompleteness.

Physical law→theory→model: Some people distinguish between laws, theories, and models, but I think they are the same thing. People who use ``law'' are likely to believe that it is objectively true and immutable; people who use ``theory'' concede that it is subject to revision; and ``model'' concedes that it is based on simplification and approximation.

Some concepts that are called ``physical laws'' are really definitions; others are, in effect, the assertion that a model predicts or explains the behavior of a system particularly well. I discuss the nature of physical models later in Think Complexity.

Determinism→indeterminism: Determinism is the view that all events are caused, inevitably, by prior events. Forms of indeterminism include randomness, probabilistic causation, and fundamental uncertainty. We come back to this topic later in the book.

These trends are not universal or complete, but the center of opinion is shifting along these axes. As evidence, consider the reaction to Thomas Kuhn's The Structure of Scientific Revolutions, which was reviled when it was published and now considered almost uncontroversial.

These trends are both cause and effect of complexity science. For example, highly abstracted models are more acceptable now because of the diminished expectation that there should be unique correct model for every system. Conversely, developments in complex systems challenge determinism and the related concept of physical law.

The excerpts so far have been from Chapter 1 of Think Complexity. Future excerpts will go into some of these topics in more depth. In the meantime, you might be interested in this timeline of complexity science (from Wikipedia):

Freshman hordes even more godless!

2012-01-29T09:38:00.000-08:00

[This is an update of an article I wrote last year, "Freshman hordes more godless than ever."]

For several years I have been following one of the most under-reported stories of the decade: the fraction of college freshmen who report no religious preference has more than tripled since 1985, from 8% to 25%, and the trend is accelerating.

Similarly, students reporting that in the last year they have never attended a religious service has grown from 8% to more than 27%.

My analysis is based on survey results from the Cooperative Institutional Research Program (CIRP) of the Higher Education Research Insitute (HERI). In 2011, more than 200,000 students at 270 colleges and universities completed the CIRP Freshman Survey, which includes questions about students’ backgrounds, activities, and attitudes.

In one question, students select their “current religious preference,” from a choice of seventeen common religions, “Other religion,” or “None.”

Another question asks students how often they “attended a religious service” in the last year. The choices are “Frequently,” “Occasionally,” and “Not at all.” Students are instructed to select “Occasionally” if they attended one or more times.

This figure shows students' responses over the history of the survey:

It's clear that both measurements are increasing, and it looks like they might be accelerating. To make that more precise, I fit a parabola to each curve. This figure shows the data for "No religion" and a least squares fit:

The red line shows the fitted model; the dark gray area shows the sample error of the model. The lighter gray area shows the sum of the sampling error and the residual error. R² for this model is 0.95; the p-values for the model and the parameters are < 0.001.

Similarly, here is the data for "No attendance" and a least squares fit:

R² for this model is 0.89, and again the p-values are effectively 0. [Note for stats geeks: as it happens, the coefficient of the linear term is close to zero, so it is not statistically significant. My first thought was to remove it from the model, but if I did that, I would understate the sampling error of the model. In this case, it is correct to keep an "insignificant" variable in the model; the fact that it is near zero doesn't mean we can ignore it, or the error associated with it.]

To test more explicitly whether the growth is accelerating, I computed the change from one year to the next in percentage points. The following figure shows these changes and a least squares fit:

Subtraction amplifies noise, so for this model R² is only 0.24, but the p-value of the slope is 0.002, so this data provides strong evidence of acceleration. The slope is 0.035 percentage points per year.

Based on this model, the predicted change for next year is 1 percentage point, so we expect the fraction of freshmen reporting no religious preference to be 26%.

The gender gap

Since the beginning of the CIRP survey, more men than women have reported no religious preference:

And the gender gap seems to be growing. Here is the difference between men and women, in percentage points, and a least squares fit:

R² for this model is 0.44; the slope is 0.035 percentage points per year, with p-value < 0.0001.

Discussion

I first wrote about this in 2007, in this article for Free Inquiry magazine. There I wrote:

College students are hardly a random sample of the population. People with more education are less likely to believe in heaven, the devil, miracles, and the literal truth of the Christian Bible. However, contrary to many people’s expectations, educated people are more likely to attend services. So, we expect the students in this sample to be less believing than the general population but also more observant.

There is reason to think that the rate of secularization in the general population is faster than what we see in this sample. Over the lifetime of the CIRP survey, college education has democratized; the percentage of high-school graduates who enter college immediately after graduation has increased from roughly 50 percent in 1970 to 65 percent in 2003. Over this time, CIRP has included more poor students, more racial minorities, and more students from families with less education. These groups tend to be more religious than the general population, so we expect their participation to increase the religiosity in the sample. Thus, the observed decrease probably underestimates the trend in the general population.

The theory of secularization—that there is a global, long-term trend away from religion—is controversial. Early sociologists, notably Max Weber, hypothesized that secularization is a predictable effect of rationalization—the increasing tendency for social actions to be based on reason rather than emotion or tradition.

In the 1950s and 1960s, many sociologists of religion defended strong theories of secularization, but, since then, several of them—including Peter Berger and Harvey Cox—have reversed their positions, arguing that religion is resurging in some areas, including the United States.

The data presented here speak directly to this debate. The CIRP survey has posed almost the same questions to a large sample of a consistently defined group for almost forty years, and the results show a clear and consistent trend away from both identification with religious sects and participation in religious services. These data make a strong case for secularization in the United States that has, if anything, accelerated in the last decade.

Data Source

Data from the 2011 CIRP Survey are reported in The American Freshman: National Norms for Fall 2011, Pryor, J. H., DeAngelo, L., Palucki Blake, L., Hurtado, S., & Tran, S., Jan 2012.

This and all previous reports are available from the HERI publications page.

Think Complexity, Part Two

2012-01-27T04:43:00.000-08:00

I repeated the simulation dozens of times, for other random initial conditions and for other numbers of oscillators. Sync every time. [...] The challenge now was to prove it. Only an ironclad proof would demonstrate, in a way that no computer ever could, that sync was inevitable; and the best kind of proof would clarify why it was inevitable.

Strogatz is a mathematician, so his enthusiasm for proofs is understandable, but his proof doesn't address what is, to me, the most interesting part the phenomenon. In order to prove that ``sync was inevitable,'' Strogatz makes several simplifying assumptions, in particular that each firefly can see all the others.

In my opinion, it is more interesting to explain how an entire valley of fireflies can synchronize despite the fact that they cannot all see each other. Think Complexity discusses how this kind of global behavior emerges from local interactions. Explanations of these phenomena often use agent-based models, which explore (in ways that would be difficult or impossible with mathematical analysis) the conditions that allow or prevent synchronization.

I am a computer scientist, so my enthusiasm for computational models is probably no surprise. I don't mean to say that Strogatz is wrong, but rather that people disagree about what questions to ask and what tools to use to answer them. These decisions are based on value judgments, so there is no reason to expect agreement. Nevertheless, there is rough consensus among scientists about which models are considered good science, and which others are fringe science, pseudoscience, or not science at all.

I claim, and this is a central thesis of the book, that the criteria this consensus is based on change over time, and that the emergence of complexity science reflects a gradual shift in these criteria.

The axes of scientific models

I have described classical models as based on physical laws, expressed in the form of equations, and solved by mathematical analysis; conversely, models of complexity systems are often based on simple rules and implemented as computations. We can think of this trend as a shift over time along two axes:

Equation-based→ simulation-based
Analysis→computation

The new kind of science is different in several other ways:

Continuous→discrete: Classical models tend to be based on continuous mathematics, like calculus; models of complex systems are often based on discrete mathematics, including graphs and cellular automata.

Linear→non-linear: Classical models are often linear, or use linear approximations to non-linear systems; complexity science is more friendly to non-linear models. One example is chaos theory, which explores the dynamics of systems of non-linear equations.

Deterministic→stochastic: Classical models are usually deterministic, which may reflect underlying philosophical determinism; complex models often feature randomness.

Abstract→detailed: In classical models, planets are point masses, planes are frictionless, and cows are spherical (see http://en.wikipedia.org/wiki/Spherical_cow). Simplifications like these are often necessary for analysis, but computational models can be more realistic.

One, two→many: In celestial mechanics, the two-body problem can be solved analytically; the three-body problem cannot. Where classical models are often limited to small numbers of interacting elements, complexity science works with larger complexes (which is where the name comes from).

Homogeneous→composite: In classical models, the elements tend to be interchangeable; complex models more often include heterogeneity.

These are generalizations, so we should not take them too seriously. And I don't mean to deprecate classical science. A more complicated model is not necessarily better; in fact, it is usually worse.

Also, I don't mean to say that these changes are abrupt or complete. Rather, there is a gradual migration in the frontier of what is considered acceptable, respectable work. Some tools that used to be regarded with suspicion are now common, and some models that were widely accepted are now regarded with scrutiny. For example, when Appel and Haken proved the four-color theorem in 1976, they used a computer to enumerate 1,936 special cases that were, in some sense, lemmas of their proof. At the time, many mathematicians did not consider the theorem truly proved. Now computer-assisted proofs are common and generally (but not universally) accepted.

Conversely, a substantial body of economic analysis is based on a model of human behavior called ``Economic man,'' or, with tongue in cheek, Homo economicus. Research based on this model was highly-regarded for several decades, especially if it involved mathematical virtuosity. More recently, this model is treated with more skepticism, and models that include imperfect information and bounded rationality are hot topics.

At this point I have laid out a lot of ideas for one article, and explained them very briefly. Think Complexity gets into these topics in more detail, but I will stop here for now. Next time I will talk about related shifts in engineering and (a little farther afield) in ways of thinking.

Think Complexity

2012-01-23T11:42:00.000-08:00

Think Complexity

On Friday I turned in the manuscript for Think Complexity, which will be published by O'Reilly Media in March (or, at least, that's the schedule). For people who can't stand to wait that long, I am going to publish excerpts here. And if you really can't wait, you can read the free version at thinkcomplex.com.

The book is about data structures and algorithms, intermediate programming in Python, computational modeling and the philosophy of science:

Data structures and algorithms: A data structure is a collection of data elements organized in a way that supports particular operations. For example, a Python dictionary organizes key-value pairs in a way that provides fast mapping from keys to values, but mapping from values to keys is slower. An algorithm is a mechanical process for performing a computation. Designing efficient programs often involves the co-evolution of data structures and the algorithms that use them. For example, in the first few chapters I present graphs, data structures that implement graphs, and graph algorithms based on those data structures.
Python programming: Think Complexity picks up where Think Python leaves off. I assume that you have read that book or have equivalent knowledge of Python. I try to emphasize fundamental ideas that apply to programming in many languages, but along the way I present some useful features that are specific to Python.
Computational modeling: A model is a simplified description of a system used for simulation or analysis. Computational models are designed to take advantage of cheap, fast computation.
Philosophy of science: The experiments and results in this book raise questions relevant to the philosophy of science, including the nature of scientific laws, theory choice, realism and instrumentalism, holism and reductionism, and epistemology.

The book is also about complexity science, which is an interdisciplinary field---at the intersection of mathematics, computer science and natural science---that focuses on discrete models of physical systems. In particular, it focuses on complex systems, which are systems with many interacting components.

Complex systems include networks and graphs, cellular automata, agent-based models and swarms, fractals and self-organizing systems, chaotic systems and cybernetic systems. Here is a framework for these topics (from http://en.wikipedia.org/wiki/Complex_systems).

A new kind of science?

In 2002 Stephen Wolfram published A New Kind of Science where he presents his and others' work on cellular automata and describes a scientific approach to the study of computational systems. There's a chapter about cellular automata in the book, but for now I'm going to borrow his title for something a little broader.

I think complexity is a ``new kind of science'' not because it applies the tools of science to a new subject, but because it uses different tools, allows different kinds of work, and ultimately changes what we mean by ``science.''

To demonstrate the difference, I'll start with an example of classical science: suppose someone asked you why planetary orbits are elliptical. You might invoke Newton's law of universal gravitation and use it to write a differential equation that describes planetary motion. Then you could solve the differential equation and show that the solution is an ellipse. Voila!

Most people find this kind of explanation satisfying. It includes a mathematical derivation---so it has some of the rigor of a proof---and it explains a specific observation, elliptical orbits, by appealing to a general principle, gravitation.

Let me contrast that with a different kind of explanation. Suppose you move to a city like Detroit that is racially segregated, and you want to know why it's like that. If you do some research, you might find a paper by Thomas Schelling called ``Dynamic Models of Segregation,'' which proposes a simple model of racial segregation (you can download the paper here).

Here is an outline of the paper's primary experiment:

The Schelling model of the city is an array of cells where each cell represents a house. The houses are occupied by two kinds of agents, labeled red and blue, in roughly equal numbers. About 10% of the houses are empty.
At any point in time, an agent might be happy or unhappy, depending on the other agents in the neighborhood. In one version of the model, agents are happy if they have at least two neighbors like themselves, and unhappy if they have one or zero.
The simulation proceeds by choosing an agent at random and checking to see whether it is happy. If so, nothing happens; if not, the agent chooses one of the unoccupied cells at random and moves.

If you start with a simulated city that is entirely unsegregated and run the model for a short time, clusters of similar agents appear. As time passes, the clusters grow and coalesce until there are a small number of large clusters and most agents live in homogeneous neighborhoods.

The degree of segregation in the model is surprising, and it suggests an explanation of segregation in real cities. Maybe Detroit is segregated because people prefer not to be greatly outnumbered and will move if the composition of their neighborhoods makes them unhappy.

Is this explanation satisfying in the same way as the explanation of planetary motion? Most people would say not, but why?

Most obviously, the Schelling model is highly abstract, which is to say not realistic. It is tempting to say that people are more complex than planets, but when you think about it, planets are just as complex as people (especially the ones that have people).

Both systems are complex, and both models are based on simplifications; for example, in the model of planetary motion we include forces between the planet and its sun, and ignore interactions between planets.

The important difference is that, for planetary motion, we can defend the model by showing that the forces we ignore are smaller than the ones we include. And we can extend the model to include other interactions and show that the effect is small. For Schelling's model it is harder to justify the simplifications.

To make matters worse, Schelling's model doesn't appeal to any physical laws, and it uses only simple computation, not mathematical derivation. Models like Schelling's don't look like classical science, and many people find them less compelling, at least at first. But as I will try to demonstrate, these models do useful work, including prediction, explanation, and design. One of the goals of Think Complexity is to explain how.

In the next excerpt, I will present some of the ways I think complexity science differs from classical science. But in the meantime I welcome comments on Part One.

Some of my best friends are crackpots

2012-01-20T06:13:00.000-08:00

I have a soft spot for crank science. Recently I visited Norumbega Tower, which is an enduring monument to the crackpot theories of Eben Norton Horsford, inventor of double-acting baking powder and faux history. But that's not what this article is about.

This article is about the Variability Hypothesis, which (quoth Wikipedia):

"originated in the early nineteenth century with Johann Meckel, who argued that males have a greater range of ability than females, especially in intelligence. In other words, he believed that most geniuses and most mentally retarded people are men. Because he considered males to be the 'superior animal,' Meckel concluded that females' lack of variation was a sign of inferiority."

I particularly like the last part, because I suspect that if it turns out that women are actually more variable, Meckel would take that as a sign of inferiority, too. Anyway, you will not be surprised to hear that evidence for the Variability Hypothesis is mixed at best.

Nevertheless, it came up in my class recently when we looked at data from the CDC's Behavioral Risk Factor Surveillance System (BRFSS), specifically the self-reported height of adult American men and women. The dataset includes responses from 154407 men and 254722 women. Here's what we found:

The average height for men is 178 cm; the average height for women is 163 cm. So men are taller, on average. No surprises so far.
For men the standard deviation is 7.7 cm; for women is it 7.3 cm. So in absolute terms, men's heights are more variable.
But to compare variability between groups, it is more meaningful to use the coefficient of variation (CV), which is the standard deviation divided by the mean. It is a dimensionless measure of variability relative to scale. For men CV is 0.0433; for women it is 0.0444.

That's very close, so we could conclude that this dataset provides no support for the Variability Hypothesis. But in keeping with the theme of this blog, I won't be happy until I have overthought the question. And along the way I want to use this problem to demonstrate Bayesian estimation in 2 dimensions.

Bayesian estimation in 2-D

In previous articles I have presented Bayesian estimation for the number of species in a biological sample and the age of a renal tumor. But those are both 1-D problems. Estimating the mean and variance of a sample is a 2-D problem. There are analytic solutions for several forms of the prior distribution, but as usual I will demonstrate a computational approach.

I'll present the code for the update first, and get back to the prior. Here's what the update looks like:

def LogUpdate(suite, evidence):
for hypo in suite.Values():
likelihood = LogLikelihood(evidence, hypo)
suite.Incr(hypo, likelihood)

suite is a Pmf object that maps from each hypothesis to its probability. While we're performing the update, the probabilities are unnormalized, which is why they are called likelihoods. And the whole thing happens under a log transform, because otherwise the likelihoods get so small they get rounded off to zero (see underflow).

Here's the function that computes likelihoods:

def LogLikelihood(evidence, hypo):
t = evidence
mu, sigma = hypo

total = Summation(t, mu)
return -len(t) * math.log(sigma) - total / 2 / sigma**2

evidence is the tuple of heights in centimeters; hypo is a tuple of hypothetical values for mu and sigma. Summation computes the sum of the squared deviations from mu; the return value is the log likelihood of getting this sample from a Gaussian distribution with the hypothetical parameters.

I pulled Summation out into a separate function because it gets called repeatedly with the same parameters, so I can improve the run time (a lot) by caching previous results:

def Summation(t, mu, cache={}):
try:
return cache[t, mu]
except KeyError:
total = sum((x-mu)**2 for x in t)
cache[t, mu] = total
return total

That's it for the update. Now we need a prior. One option is a uniform prior where mu can be any value and sigma can be any positive value. In practice, there is no point in evaluating values far from the true parameters, because their likelihoods will be vanishingly small. So we'll use the data to choose the range for the prior, then use the data to update the prior.

That might sound completely bogus, because we are using the same data twice. But choosing the range for the prior is just a computational convenience; it has no effect on the result, as long as the range is wide enough that the likelihood for values outside the range is effectively zero.

To choose the prior range for mu and sigma, I use the conventional estimators and their standard errors. Here's the code:

def MakeUniformPrior(t, num_points, label, spread=3.0):

# estimate mean and stddev of t
n = len(t)
xbar, S2 = thinkstats.MeanVar(t)
sighat = math.sqrt(S2)

# compute standard error for mu and the range of ms
stderr_xbar = sighat / math.sqrt(n)
mspread = spread * stderr_xbar
ms = numpy.linspace(xbar-mspread, xbar+mspread, num_points)

# compute standard error for sigma and the range of ss
stderr_sighat = sighat / math.sqrt(2 * (n-1))
sspread = spread * stderr_sighat
ss = numpy.linspace(sighat-sspread, sighat+sspread, num_points)

# populate the PMF
pmf = Pmf.Pmf(name=label)
for m in ms:
for s in ss:
pmf.Set((m, s), 1)
return ms, ss, pmf

Finally, here's the code that makes the prior, updates it, and normalizes the posterior:

def EstimateParameters(t, label, num_points=31):

xs, ys, suite = MakeUniformPrior(t, num_points, label)

suite.Log()

LogUpdate(suite, tuple(t))

suite.Exp()

suite.Normalize()

return xs, ys, suite

Ok, that was more code than I usually put in these articles. You can download it all from thinkcomplex.com/bayes_height.py.

Results

The result is a map from pairs of mu and sigma to probabilities; one simple way to display the map is a contour plot. Here is the result for men:

And for women:

As expected, the central location is near the estimated parameters. But with this result in numerical form is that it is easy to compute the marginal distributions for mu and sigma, or their confidence intervals. We can also compute posterior distributions for CV, and compare the distributions CV for men and women:

Since there is no overlap between the distributions, we can conclude that the small observed difference is statistically significant; or, if we want to be more Bayesian about it, we can compute the probability that the CV for women is higher, which is pretty close to 1.

Finally we conclude:

The standard deviation of height for men is higher.
But the coefficient of variation is slightly lower.
Even though the difference is tiny, it is very unlikely to be due to chance.
And therefore women are inferior.

Hey, I told you I have a soft spot for crank science.

Frank is a scoundrel, probably

2012-01-05T09:44:00.000-08:00

My friend Ted Bunn wrote an article, "Who knows what evil lurks in the hearts of men? The Bayesian doesn’t care," inspired in part by one of my posts, "Repeated tests: how bad can it be?"

He presents this scenario:

Frank and Betsy are wondering whether a particular coin is a fair coin (i.e., comes up heads and tails equally often when flipped). Frank, being a go-getter type, offers to do some tests to find out. He takes the coin away, flips it a bunch of times, and eventually comes back to Betsy to report his results.

“I flipped the coin 3022 times,” he says, “and it came up heads 1583 times. That’s 72 more heads than you’d expect with a fair coin. I worked out the p-value — that is, the probability of this large an excess occurring if the coin is fair — and it’s under 1%. So we can conclude that the coin is unfair at a significance level of 1% (or ’99% confidence’ as physicists often say).”

And suggests that there are two ways Betsy can interpret this report (again, quoting Ted):

Frank is an honest man, who has followed completely orthodox (frequentist) statistical procedure. To be specific, he decided on the exact protocol for his test (including, for some reason, the decision to do 3022 trials) in advance.
Frank is a scoundrel who, for some reason, wants to reach the conclusion that the coin is unfair. He comes up with a nefarious plan: he keeps flipping the coin for as long as it takes to reach that 1% significance threshold, and then he stops and reports his results.

Finally, Ted asks, "What should Betsy conclude on the basis of the information Frank has given her?" If you want to know Ted's answer, you have to read his article (and you should -- it is very interesting).

But I'm going to give you my answer: Frank is a scoundrel. Well, probably.

I'll follow Ted by making one more assumption: "Suppose that Betsy’s initial belief is that 95% of coins are fair — that is, the probability P that they come up heads is exactly 0.5." Now we can evaluate the evidence that Frank is a scoundrel.

If Frank is a scoundrel, the probability that he reports a positive result is 1, provided that he is willing to keep flipping long enough, or 1-x, for small x, if we put a bound on the number of flips he is willing to do. So

P(positive test | Frank is a scoundrel) = 1-x

If Frank is honest, then the probability of a positive result is

P(fair coin) P(false positive | fair coin) + P(biased coin) * P(true positive | biased coin)

Betsy believes that P(fair coin) is 95% and P(biased coin) is 5%. Since Frank's significance level is 1%, P(false positive | fair coin) is 1%.

The probability of a true positive is the power of the test, which depends on how biased the coin actually is. But I will make the approximation that 3022 flips is enough to detect any substantial bias, so I'll take P(true positive | biased coin) = 1-y, for small y. So,

P(positive test | Frank is honest) = (0.95)(0.01) + (0.05)(1-y)

As x and y get small, the likelihood ratio is (1 / 0.0595), which is about 17. So that is fairly strong evidence that Frank is a scoundrel.

I don't know about Betsy's prior beliefs about Frank, so you will have to fill in your own punchline about her posterior.

Are first babies more likely to be light?

2011-11-30T07:29:00.000-08:00

Way back in February I analyzed data from the National Survey of Family Growth (NSFG) and answered the question "Are first babies more likely to be late?" Now I am getting back to that data to look at the related question, "Are first babies more likely to be light?"

In Think Stats (Chapter 7), I showed that the mean birth weight for first babies is lower than the mean for others by about 2 ounces, and that difference is statistically significant. But there is also a relationship between birth weight and mother's age, and the mothers of first babies tend to be younger. So the question is: if we control for the age of the mother, are first babies still lighter?

Several students in my class are working on projects involving multivariate regression, so I want to use this question as an example. I will follow the steps I recommend to my students:

1) Before looking at relationships between variables, look at each variable in isolation. In particular, characterize the distributions and identify issues like outliers or long tails.

2) Look at the variables pairwise. For each pair, look at CDFs and scatterplots, and compute correlations and/or least squares fits.

3) If there seem to be relationships among the variables, look for ways to separate the effects, either by breaking the data into subsets or doing multivariate regression.

So let's proceed.

One variable at a time

The variables I'll use are

Mother's age in years,
Birth weight in ounces, and
first, which is a dummy variable, 1 for first babies and 0 for others.

For live births we have 9148 records with valid ages; here is the distribution:

It looks like there is a little skew to the right, but other than that, nothing to worry about.

We have 9038 records with valid weights, with this distribution:

The middle of the distribution is approximately normal, with some jaggies due to round-off. In the tails there are some values that are are certainly errors, but it is hard to draw a clear line between exceptional cases and bad data. It might be a good idea to exclude some extreme values, but for the analysis below I did not.

There are only two values for first, so it's not much of a distribution, but with categorical data, we should check that we have enough values in each bin. As it turns out, there are 4413 first babies and 4735 others, so that's just fine.

One pair at a time

To characterize the relationship between weight and first, we compare the CDF of weights for first babies and others:

There is some space between the distributions, and the difference in means is 2 ounces. It's not obvious whether that difference is statistically significant, but it is, with p < 0.001.

Similarly we can compare the CDF of mother's age for first babies and others:

Here there is clearly space between the distributions. The difference in means is 3.6 years, which is significant (no surprise this time).

It is not as easy to see the relationship between age and birthweight. I often tell my students to start with a scatterplot:

But in this case it's not much help. If there is a relationship there, it is hard to see. An alternative is to break the age range into bins and compute the mean in each bin. Here's what that looks like with 2-year bins:

We can't take the first and last points too seriously; there are not many cases in those bins. And ideally I should represent the variability in each bin so we have a sense of whether the apparent differences are real. Based on this figure, it looks like there is a relationship, but it might be nonlinear.

Pearson's coefficient of correlation between age and birthweight is 0.07, which is small but statistically significant. Spearman's coefficient of correlation is 0.10; the difference between the two coefficients is another warning that the relationship is nonlinear.

If we ignore the non-linearity for now, we can compute a least squares fit for birthweight as a function of age. The slope is 0.28, which means that we expect an additional 0.28 ounces per year of age. Since first mothers are 3.6 years younger than others, we expect their babies to be 1.0 ounces lighter. In fact, they are 2.0 ounces lighter, so the linear model of weight vs age accounts for 50% of the observed difference.

Multiple regression

So far I have been using the thinkstats libraries to compute correlation coefficients and least squares fits. But for multiple regression I am going to break out rpy, which is the Python interface to R, a "language and environment for statistical computing and graphics. It is a GNU project which is similar to the S language and environment which was developed at Bell Laboratories." To see how it works, you can download the code I used in this section: age_lm.py

The first model I ran is the same linear model we were just looking at:

weight = intercept + slope * age

In R's shorthand, that's:

weights ~ ages

Here are the results:

Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.28635 1.08775 100.470 < 2e-16 ***
ages 0.27926 0.04258 6.559 5.72e-11 ***

Multiple R-squared: 0.004738, Adjusted R-squared: 0.004628

The intercept is 109 ounces; the slope is 0.28 ounces per year, which what we saw before---comparing my implementation to R makes me more confident that R is correct :)

The standard error provides a confidence interval for the estimates; plus or minus 2 standard errors is (roughly) a 95% confidence interval.

The last column is the p-value, which is very small, indicating that the slope and intercept are significantly different from 0. The t-value is the test statistic used to compute the p-value, but I don't know why it gets reported; it doesn't mean much.

Since the p-values are so small, it might be surprising that R² is so low, only 0.0047. But there is no contradiction. We can say with high confidence that there is a relationship between these variables; nevertheless, birth weight is highly variable, and even if you know the age of the mother, that does not reduce the variability by much.

To understand adjusted R², consider this: if you add more explanatory variables to a model, R² usually goes up even if there is no real relationship. The adjusted R² takes this into account, which makes it more meaningful to compare models with a different number of variables.

Now let's add first as an explanatory variable, so the model looks like this:

weights ~ first + ages

Here are the results:

Estimate Std. Error t value Pr(>|t|)
(Intercept) 110.62791 1.24198 89.073 < 2e-16 ***
first -1.11688 0.49940 -2.236 0.0253 *
ages 0.24708 0.04493 5.499 3.93e-08 ***

Multiple R-squared: 0.005289, Adjusted R-squared: 0.005069

The coefficient for first is -1.1, which means that we expect first babies to be 1.1 ounces lighter, controlling for age. The p-value for this estimate is 2.5%, which I consider borderline significant.

The coefficient for ages is about the same as before, and significant. And R² is a little higher, but still small.

But remember that the relationship between weight and age is non-linear. We can explore that by introducing a new variable:

ages2 = ages^2

Now if we run this model:

weights ~ ages + ages2

We are effectively fitting a parabola to the weight vs age curve.

Estimate Std. Error t value Pr(>|t|)
(Intercept) 89.151237 4.407677 20.226 < 2e-16 ***
ages 1.898151 0.346071 5.485 4.25e-08 ***
ages2 -0.031002 0.006577 -4.714 2.47e-06 ***

Multiple R-squared: 0.00718, Adjusted R-squared: 0.00696

Estimates for both variables are significant. The coefficient for ages2 is negative, which means that the parabola has downward curvature, as expected. And R² is a little bigger.

Finally, we can bring it all together:

weights ~ first + ages + ages2

With this model we get:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 91.07627 4.56813 19.937 < 2e-16 ***

first -0.80708 0.50373 -1.602 0.109

ages 1.79807 0.35163 5.113 3.23e-07 ***

ages2 -0.02953 0.00664 -4.447 8.80e-06 ***

Multiple R-squared: 0.007462, Adjusted R-squared: 0.007132

When we include the parabolic model of weight and age, the coefficient for first gets smaller, and the p-value is 11%.

I conclude that the difference in weight for first babies is explained by the difference in mothers' ages. When we control for age, the difference between first babies and others is no longer statistically significant. It is possible that there is a small difference in weight for first babies, but this dataset provides little evidence for it.

Estimating the age of renal tumors

2011-11-28T05:11:00.000-08:00

A few weeks ago I read this post on reddit.com/r/statistics:

"I have Stage IV Kidney Cancer and am trying to determine if the cancer formed before I retired from the military. ... Given the dates of retirement and detection is it possible to determine when there was a 50/50 chance that I developed the disease? Is it possible to determine the probability on the retirement date? My tumor was 15.5 cm x 15 cm at detection. Grade II."

I contacted the original poster and got more information; I learned that veterans get different benefits if it is "more likely than not" that a tumor formed while they were in military service (among other considerations).

Because renal tumors grow slowly, and often do not cause symptoms, they are often left untreated. As a result, we can observe the rate of growth for untreated tumors by comparing scans from the same patient at different times. Several papers have reported these growth rates.

I collected data from a paper by Zhang et al. I contacted the authors to see if I could get raw data, but they refused on grounds of medical privacy. Nevertheless, I was able to extract the data I needed by printing one of their graphs and measuring it with a ruler. It's silly, but it works.

They report growth rates in reciprocal doubling time (RDT), which is in units of doublings per year. So a tumor with RDT=1 doubles in volume each year; with RDT=2 it quadruples in the same time, and with RDT=-1, it halves. The following figure shows the distribution of RDT for 53 patients:

The squares are the data points from the paper; the line is a model I fit to the data. The positive tail fits an exponential distribution well, so I used a mixture of two exponentials.

As a simple model of tumor growth, I chose the median value of RDT, which is 0.45, and used that to estimate the age of a tumor with maximum dimension 15.5 cm. Here's what I wrote in my letter to the Veterans Benefits Administration:

In the largest study I reviewed (53 patients) the median volume doubling time is 811 days. By definition of median, 50% of observed tumors grew faster and 50% slower. By geometry, the doubling time for the maximum linear dimension is approximately (811)(3) = 2433 days or 6.7 years. Therefore, for a tumor with maximum linear dimension 15.5 cm on [diagnosis date], it is as likely as not that the size on [discharge date] was 6 cm.
If the diameter of the tumor on [discharge date] were 1 mm and it grew to 15.5 cm by [diagnosis date], the effective volume doubling time would be 150 days. Fewer than half of the tumors in the studies I reviewed grew at this rate or faster, so it is more likely than not that the tumor grew more slowly.

Based on this analysis, I conclude that it is more likely than not that this tumor formed prior to [discharge date].

I think this model is sufficient to answer the question as posed, but it occurred to me later (in the shower, where all good ideas come from) that we can do better. By sampling from the distribution of growth rates and generating simulated tumor histories, we can estimate the distribution of size as a function of time and then, using Bayes's Theorem, get the distribution of age as a function of size.

Here's how. The simulation starts with a small tumor (0.3 cm) and runs these steps:

Choose a growth rate from the distribution of RDT.
Compute the size of the tumor at the end of an 8 month interval (that's the median interval between measurements in the data source).
Repeat until the tumor is 20 cm in diameter.

This figure shows 100 simulated growth trajectories:

The line at 10 cm shows the range of ages for tumors at that size: the fastest-growing tumor gets there in 8 years; the slowest takes more than 35.

By drawing the line at different sizes, we can estimate the distribution of age as a function of size. There's an implicit use of Bayes's Theorem in there, but because I did everything discretely, I didn't have to think too hard. This figure shows the distribution of age for a few different sizes:

Not surprisingly, bigger tumors are likely to be older. For any size, we can generate the CDF and compute the median, interquartile range, and 90% confidence interval. Here's what that looks like (with size on a log scale):

The points are data from simulation, which produces some variability due to discrete approximation. The lines are fitted to the data.

With these results, doctors can look up the size of a tumor and get the distribution of ages; for example, the median age of a 15 cm tumor is 27 years, with interquartile range 22-31 and 90% confidence interval 16-39 years.

This model yields more detail than the simple model I started with, but the results are qualitatively similar; a tumor this size is more likely than not to have formed prior to the original poster's date of discharge. It looks like there is also a good chance that it formed prior to enlistment, but I don't know what the VBA makes of that.

-----

I think this model makes the best use of the available data, but there are several limitations:

1) The factors that limit tumor growth are different for very small tumors, so the observed data doesn't apply. We can extrapolate back to when the tumor was small (I chose 0.3 cm, a bit smaller than the smallest tumor in the study). That gives us a lower bound on the age of the tumor, but we can't say much about when the first cancer cell appeared.

2) The distribution of growth rates is based on a sample of 53 patients. A different sample would yield a different distribution. I could use resampling to characterize this source of error, but haven't.

3) The growth model does not take into account tumor subtype or grade, which is consistent with the conclusion of Zhang et al: “Growth rates in renal tumors of different sizes, subtypes and grades represent a wide range and overlap substantially.”

4) In our model of tumor growth, the growth rate during each interval is independent of previous growth rates. It is plausible that, in reality, tumors that have grown quickly in the past are more likely to grow quickly.

If this correlation exists, it affects the location and spread of the results. For example, running simulations with ρ = 0.4 increases the estimated median age by about a year, and the interquartile range
by about 3 years. However, if there were a strong serial correlation in growth rate, there would be also be a correlation between tumor volume and growth rate, and prior work has shown no such relationship.

There could still be a weak serial correlation, but since there is currently no evidence for it, I ran these simulations with ρ = 0.

Comment on "Racism and Meritocracy"

2011-11-21T11:39:00.000-08:00

WARNING: This article is on a topic that elicits emotional reactions. I welcome comments, but please make them thoughtful and keep them civil.

Eric Ries wrote an article for TechCrunch last week, talking about racism and meritocracy among Silicon Valley entrepreneurs. It's a good article; you should read it and then come back.

Although I mostly agree with him, Ries undermines his argument with a statistical bait-and-switch: he starts out talking about race, but most of the article (and the slide deck he refers to) are about gender. Unfortunately, for both his argument and the world, the race gap is bigger than the gender gap, and it is compounded because racial minorities, unlike women, are minorities.

To quantify the size of the gap, I use data from Academically Adrift, a recent book that reports the results from the Collegiate Learning Assessment (CLA) database, collected by the Council for Aid to Education, "a national nonprofit organization ... established in 1952 to advance corporate support of education and to conduct policy research on higher education..."

Here is a description of the CLA:

"The CLA consists of three types of prompts within two types of task: the Performance Task and the Analytic Writing Task...The Analytic Writing Task includes a pair of prompts called Make-an-Argument and Critique-an-Argument.

"The CLA uses direct measures of skills in which students perform cognitively demanding tasks... All CLA measures are administered online and contain open-ended prompts that require constructed responses. There are no multiple-choice questions. The CLA tasks require that students integrate critical thinking and written communication skills. The holistic integration of these skills on the CLA tasks mirrors the requirements of serious thinking and writing tasks faced in life outside of the classroom. "

This is not your father's SAT. The exam simulates realistic workplace tasks and assesses skills that are relevant to many jobs, including (maybe especially) entrepreneurship.

On this assessment, the measured differences between black and white college students are stark. For white college students, the mean and standard deviation are 1170 ± 179. For black students, they are 995 ± 167.

To get a sense of what that difference looks like, suppose there are just two groups, which I call "blue" and "green" as a reminder that I am presenting an abstract model and not a realistic description. This figure shows Gaussian distributions with the parameters reported in Academically Adrift:

The difference in means is 175 points, which is about one standard deviation. If we select people from the upper tail, the majority are blue. But the situation is even worse if greens are a minority. If greens make up 20% of the population, the picture looks like this:

The fraction of greens in the upper tail is even smaller. If, as Ries suggests, "Here in Silicon Valley, we’re looking for the absolute best and brightest, the people far out on the tail end of aptitude," the number of greens in that tail is very small.

How small? That depends on where we draw the line. If we select people who score above 1200, which includes 37% of the population, we get 6% greens (remember that they are 20% of the hypothetical population). Above 1300 the proportion of greens is 3%, and above 1400 only 2%.

And that's not very "far out on the tail end of aptitude." Above 1500, we are still talking about 3% of the general population, but more than 99% of them are blue. So in this hypothetical world of blues and greens, perfect meritocracy does not lead to proportional representation.

Ries suggests that blind screening of applicants might help. I think the system he proposes is a good idea, because it improves fairness and also the perception of fairness. But if the racial gap in Y Combinator's applicant pool is similar to the racial gap in CLA scores, making the selection process more meritocratic won't make a big difference.

These numbers are bad. I'm sorry to be reporting them, and if I know the Internet, some people are going to call me a racist for doing it. But I didn't make them up, and I'm pretty sure I did the math right. Of course, you are welcome to disagree with my conclusions.

Here are some of the objections I expect:

1) The CLA does not capture the full range of skills successful entrepreneurs need.

Of course it doesn't; no test could. But I chose the CLA because I think it assesses thinking skills better than other standardized tests, and because the database includes "over 200,000 student results across hundreds of colleges." I can't think of a better way to estimate the magnitude of the racial gap in the applicant pool.

2) The application process is biased against racial minorities and women.

The statistics I am reporting here, and my analysis of them, don't say anything about whether or not the application process is biased. But they do suggest (a) We should not assume that because racial minorities are underrepresented among Silicon Valley entrepreneurs, racial bias explains a large part of the effect, and (b) We should not assume that eliminating bias from the process will have a large effect.

Of course, trying to eliminate bias is the right thing to do, whether the effect is big or small.

-----

NOTE: The range of scores for the CLA was capped at 1600 until 2007, which changed the shape of the distribution at the high end. For those years, the Gaussian distributions in the figures are not exactly right, but I don't think it affects my analysis much. Since 2007, scores are no longer capped, but I don't know what the tail of the distribution looks like now.

EDIT 11-28-11: I revised a few sentences to clarify whether I was talking about representation or absolute numbers. The fraction of greens in the population affects the absolute numbers in the tail but not their representation.

Girl Named Florida solutions

2011-11-10T11:09:00.000-08:00

In The Drunkard's Walk, Leonard Mlodinow presents "The Girl Named Florida Problem":

"In a family with two children, what are the chances, if one of the children is a girl named Florida, that both children are girls?"

I like this problem, and I use it on the first day of my class to introduce the topic of conditional probability. But I've decided that it's too easy. To give it a little more punch, I've decided to combine it with the Red-Haired Problem from last week:

In a family with two children, what are the chances, if at least one of the children is a girl with red hair, that both children are girls?

Just like last week, you can make some simplifying assumptions:

About 2% of the world population has red hair. You can assume that the alleles for red hair are purely recessive. Also, you can assume that the Red Hair Extinction theory is false, so you can apply the Hardy–Weinberg principle. And you can ignore the effect of identical twins.

Before I present my solution, I want to sneak up on it with a series of warm-up problems.

P[GG | two children]: if a family has two children, what is the chance that they have two girls?
P[GG | two children, at least one girl]: if we know they have at least one girl, what is the chance that they have two girls?
P[GG | two children, older child is a girl]: if the older child is a girl, what is the chance that they have two girls?
P[GG | two children, at least one is a girl named Florida].
P[GG | two children, at least one is a girl with red hair, and the parents have brown hair].
P[GG | two children, at least one is a girl with red hair].

Problem 1: P[GG | two children]

If we assume that the probability that each child is a girl is 50%, then P[GG | two children] = 1/4.

Problem 2: P[GG | two children, at least one girl]

There are four equally-likely kinds of two child families: BB, BG, GB and GG. We know that BB is out, so the conditional probability is

P[GG | at least one girl] = P[GG and at least one girl] / P[at least one girl] = 1/3.

Problem 3: P[GG | two children, older child is a girl]

Now there are only two possible families, GB and GG, so the conditional probability is 1/2. Informally we can argue that once we know about the older child we can treat the younger child as independent. But if there's one thing we learn from this problem, it's that our intuition for independence is not reliable.

Problem 4: P[GG | two children, at least one girl named Florida]

Here's the one that makes people's head hurt. For each child, there are three possibilities, boy, girl not named Florida, and girl named Florida, with these probabilities:

B: 1/2
G: 1/2 - x
GF: x

where x is the unknown percentage of people who are girls named Florida. Of families with at least one girl named Florida, there are these possible combinations, with these probabilities

B GF: 1/2 x
GF B: 1/2 x
G GF: x (1/2 - x)
GF G: x (1/2 - x)
GF GF: x^2

The highlighted cases have two girls, so the probability we want is the sum of the highlighted cases over the sum of all cases. With a little algebra, we get:

P(GG | at least one girl named Florida) = (1 - x) / (2 - x)

Assuming that Florida is not a common name, x approaches 0 and the answer approaches 1/2. So it turns out, surprisingly, that the name of the girl is relevant information.

As x approaches 1/2, the answer converges on 1/3. For example, if we know that at least one child is a girl with two X chromosomes, x is close to 1/2 and the problem reduces to Problem 2.

If this problem is still making your head hurt, this figure might help:

Here B a boy, Gx is a girl with some property X, and G is a girl who doesn't have that property. If we select all families with at least one Gx, we get the five blue squares (light and dark). The families with two girls are the three dark blue squares.

If property X is common, the ratio of dark blue to all blue approaches 1/3. If X is rare, the same ratio approaches 1/2.

Problem 5: P[GG | two children, at least one girl with red hair, parents have brown hair]

If the parents have brown hair and one of their children has red hair, we know that both parents are heterozygous, so their chance of having a red-haired girl is 1/8.

Using the girl-named-Florida formula, we get

P[GG | two children, at least one girl with red hair, parents have brown hair] = (1 - 1/8) / (2 - 1/8) = 7/15.

And finally:

Problem 6: P[GG | two children, at least one girl with red hair]

In this case we don't know the genotype of the parents. There are three possibilities: Aa Aa, Aa aa, and aa aa.

We follow these steps:

1) Use the prevalence of red hair to compute the prior probabilities of each parental genotype.

2) Use the evidence (at least one girl with red hair) to compute the posteriors.

3) For each combination, compute the conditional probability. We have already computed

P(GG | two children, at least one with red hair, Aa Aa) = 7/15

The others are

P(GG | two children, at least one with red hair, Aa aa) = 3/7
P(GG | two children, at least one with red hair, aa aa) = 1/3

4) Apply the law of total probability to get the answer.

I'm too lazy to do the algebra, so I got Mathematica to do it for me. Here is the notebook with the answer.

In general, if p is the prevalence of red hair alleles,

P(GG | two children, at least one with red hair) = (p^2 + 2p - 7) / (p^2 + 2p - 15)

If the prevalence of red hair is 0.02, then p = sqrt(0.02) = 0.141, and

P(GG | two children, at least one with red hair) = 45.6%

Congratulations to Professor Ted Bunn at the University of Richmond, the only person who submitted a correct answer before the deadline!

At least, I think it's the right answer. Maybe we both made the same mistake.

For more fun with probability, see Chapter 5 of my book, Think Stats, which you can read here, or buy here.

The red-haired girl named Florida

2011-11-07T13:12:00.000-08:00

In The Drunkard's Walk, Leonard Mlodinow presents "The Girl Named Florida Problem":

"In a family with two children, what are the chances, if one of the children is a girl named Florida, that both children are girls?"

In a family with two children, what are the chances, if at least one of the children is a girl with red hair, that both children are girls?

[Edit: I clarified the wording to say that at least one of the children is a girl with red hair. To be even more precise, I am looking for the conditional probability P(two girls | at least one girl with red hair).]

Just like last week, you can make some simplifying assumptions:

About 2% of the world population has red hair. You can assume that the alleles for red hair are purely recessive. Also, you can assume that the Red Hair Extinction theory is false, so you can apply the Hardy–Weinberg principle. And you can ignore the effect of identical twins.

Here is my solution.

For more fun with probability, see Chapter 5 of my book, Think Stats, which you can read here, or buy here.

Somebody bet on the Bayes

2011-11-03T08:55:00.000-07:00

In last week's post I wrote solutions to some of my favorite Bayes's Theorem problems, and posed this new problem:

If you meet a man with (naturally) red hair, what is the probability that neither of his parents has red hair?

Hints: About 2% of the world population has red hair. You can assume that the alleles for red hair are purely recessive. Also, you can assume that the Red Hair Extinction theory is false, so you can apply the Hardy–Weinberg principle.

Given the prevalence of red hair, we know what fraction of the population is homozygous recessive. To solve the problem, we also need to know how many are heterozygous and homozygous dominant.

I'll use a to represent the recessive allele (or alleles) for red hair, and A for the dominant alleles that code for other colors. p is the prevalence of a and q is the prevalence of A, so p+q = 1.

If these prevalences are not changing over time, and they don't affect people's mating decisions, we can invoke the Hardy-Weinberg principle to get:

P(AA) = prevalence of homozygous dominant = q**2
P(Aa) = prevalence of heterozygous = 2 * p * q
P(aa) = prevalence of homozygous recessive = p**2

And P(AA) + P(Aa) + P(aa) = 1.

Given (aa), we can compute p, q and:

P(Aa) = 0.243
P(AA) = 0.737

Now if a child has red hair, both parents have at least one recessive allele, so the possible combinations are (Aa Aa), (Aa aa), and (aa aa). If we assume, again, that mating decisions are not based on hair color, we can get the prevalence of each parental pair:

P(Aa Aa) = Aa**2

P(Aa aa) = 2 Aa aa

P(aa aa) = aa**2

And we can use those as the priors. The evidence is

E: the child has red hair

To get the likelihoods, we apply Mendelian genetics:

P(E | Aa Aa) = 0.25

P(E | Aa aa) = 0.5

P(E | aa aa) = 1.0

Finally, applying Bayes's Theorem, we have

P(Aa Aa | E) = P(Aa Aa) P(E | Aa Aa) / P(E)

And the answer is 0.737. With a little algebra we can show that

P(Aa Aa | E) = P(AA)

That is, the probability that neither parent has red hair is exactly the fraction of the population that is homozygous dominant.

Almost 75% of red-haired people come from parents with non-red hair, which might explain why a "red haired child" is a metaphor for a child that doesn't resemble his parents. In the expression, "beat like a red haired step child," some of the humor (for people who find child abuse funny) comes from the suggestion that the parentage of a red haired child is suspect.

But why should red hair be funnier than blue eyes or blonde hair? It turns out that we can answer this question mathematically. If the prevalence of red hair were higher, say 10%, most red haired people would have at least one red-haired parent, and that would be less funny.

In general, as the prevalence of the recessive phenotype increases, the potential for amusing insinuations of infidelity decreases; near 0 it drops off steeply, as shown in this figure:

Since I believe I am the first person to quantify this effect, I humbly submit that it should be called "Downey's inverse law of mailman jokes."

For more fun with probability, see Chapter 5 of my book, Think Stats, which you can read here, or buy here.

-----

If you don't get the title of this post, it is a play on "Somebody bet on the bay," a lyric from the minstrel song "Camptown Races." A bay is a horse with a reddish-brown coat, so I thought it was a pretty good fit.

All your Bayes are belong to us!

2011-10-27T07:45:00.000-07:00

This week's post contains solutions to My Favorite Bayes's Theorem Problems, and one new problem. If you missed last week's post, go back and read the problems before you read the solutions!

If you don't understand the title of this post, brush up on your memes.

1) The first one is a warm-up problem. I got it from Wikipedia (but it's no longer there):

Suppose there are two full bowls of cookies. Bowl #1 has 10 chocolate chip and 30 plain cookies, while bowl #2 has 20 of each. Our friend Fred picks a bowl at random, and then picks a cookie at random. We may assume there is no reason to believe Fred treats one bowl differently from another, likewise for the cookies. The cookie turns out to be a plain one. How probable is it that Fred picked it out of Bowl #1?

First the hypotheses:
A: the cookie came from Bowl #1
B: the cookie came from Bowl #2

And the priors:
P(A) = P(B) = 1/2

The evidence:
E: the cookie is plain

And the likelihoods:
P(E|A) = prob of a plain cookie from Bowl #1 = 3/4
P(E|B) = prob of a plain cookie from Bowl #2 = 1/2

Plug in Bayes's theorem and get
P(A|E) = 3/5

You might notice that when the priors are equal they drop out of the BT equation, so you can often skip a step.

2) This one is also an urn problem, but a little trickier.

The blue M&M was introduced in 1995. Before then, the color mix in a bag of plain M&Ms was (30% Brown, 20% Yellow, 20% Red, 10% Green, 10% Orange, 10% Tan). Afterward it was (24% Blue , 20% Green, 16% Orange, 14% Yellow, 13% Red, 13% Brown).

A friend of mine has two bags of M&Ms, and he tells me that one is from 1994 and one from 1996. He won't tell me which is which, but he gives me one M&M from each bag. One is yellow and one is green. What is the probability that the yellow M&M came from the 1994 bag?

Hypotheses:
A: Bag #1 from 1994 and Bag #2 from 1996
B: Bag #2 from 1994 and Bag #1 from 1996

Again, P(A) = P(B) = 1/2.

The evidence is:
E: yellow from Bag #1, green from Bag #2

The likelihoods are
P(E|A) = (0.2)(0.2)
P(E|B) = (0.1)(0.14)

So P(A|E) = 40 / 54 ~ 0.74

By introducing the terms Bag #1 and Bag #2, rather than "the bag the yellow M&M came from" and "the bag the green came from," I avoided the part of this problem that can be tricky: keeping the hypotheses and the evidence straight.

3) This one is from one of my favorite books, David MacKay's Information Theory, Inference, and Learning Algorithms:

Elvis Presley had a twin brother who died at birth. What is the probability that Elvis was an identical twin?

To answer this one, you need some background information: According to the Wikipedia article on twins: ``Twins are estimated to be approximately 1.9% of the world population, with monozygotic twins making up 0.2% of the total---and 8% of all twins.''

There are several ways to set up this problem; I think the easiest is to think about twin birth events, rather than individual twins, and to take the fact that Elvis was a twin as background information.

So the hypotheses are
A: Elvis's birth event was an identical birth event
B: Elvis's birth event was a fraternal twin event

If identical twins are 8% of all twins, then identical birth events are 8% of all twin birth events, so the priors are

P(A) = 8%
P(B) = 92%

The relevant evidence is
E: Elvis's twin was male

So the likelihoods are
P(E|A) = 1
P(E|B) = 1/2

Because identical twins are necessarily the same sex, but fraternal twins are equally likely to be opposite sex (or, at least, I assume so). So

P(A|E) = 8/54 ~ 0.15.

The tricky part of this one is realizing that the sex of the twin provides relevant information!

4) Also from MacKay's book:

Two people have left traces of their own blood at the scene of a crime. A suspect, Oliver, is tested and found to have type O blood. The blood groups of the two traces are found to be of type O (a common type in the local population, having frequency 60%) and of type AB (a rare type, with frequency 1%). Do these data (the blood types found at the scene) give evidence in favour [sic] of the proposition that Oliver was one of the two people whose blood was found at the scene?

For this problem, we are not asked for a posterior probability; rather we are asked whether the evidence is incriminating. This depends on the likelihood ratio, but not the priors.

The hypotheses are
X: Oliver is one of the people whose blood was found
Y: Oliver is not one of the people whose blood was found

The evidence is
E: two blood samples, one O and one AB

We don't need priors, so we'll jump to the likelihoods. If X is true, then Oliver accounts for the O blood, so we just have to account for the AB sample:

P(E|X) = 0.01

If Y is true, then we assume the two samples are drawn from the general population at random. The chance of getting one O and one AB is

P(E|Y) = 2(0.6)(0.01) = 0.012

Notice that there is a factor of two here because there are two permutations that yield E.

So the evidence is slightly more likely under Y, which means that it is actually exculpatory! This problem is a nice reminder that evidence that is consistent with a hypothesis does not necessarily support the hypothesis.

5) I like this problem because it doesn't provide all of the information. You have to figure out what information is needed and go find it.

According to the CDC, ``Compared to nonsmokers, men who smoke are about 23 times more likely to develop lung cancer and women who smoke are about 13 times more likely.''
If you learn that a woman has been diagnosed with lung cancer, and you know nothing else about her, what is the probability that she is a smoker?

I find it helpful to draw a tree:

If y is the fraction of women who smoke, and x is the fraction of nonsmokers who get lung cancer, the number of smokers who get cancer is proportional to 13xy, and the number of nonsmokers who get lung cancer is proportional to x(1-y).

Of all women who get lung cancer, the fraction who smoke is 13xy / (13xy + x(1-y)).

The x's cancel, so it turns out that we don't actually need to know the absolute risk of lung cancer, just the relative risk. But we do need to know y, the fraction of women who smoke. According to the CDC, y was 17.9% in 2009. So we just have to compute

13y / (13y + 1-y) ~ 74%

This is higher than many people guess.

6) Next, a mandatory Monty Hall Problem. First, here's the general description of the scenario, from Wikipedia:

Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say Door A [but the door is not opened], and the host, who knows what's behind the doors, opens Door B, which has a goat. He then says to you, "Do you want to pick Door C?" Is it to your advantage to switch your choice?

The answer depends on the behavior of the host when the car is behind Door A. In this case the host can open either B or C. Suppose he chooses B with probability p and C otherwise. What is the probability that the car is behind Door A (as a function of p)?

The hypotheses are
A: the car is behind Door A
B: the car is behind Door B
C: the car is behind Door C

And the priors are
P(A) = P(B) = P(C) = 1/3

The likelihoods are
P(E|A) = p, because in this case Monty has a choice and chooses B with probability p,
P(E|B) = 0, because if the car were behind B, Monty would not have opened B, and
P(E|C) = 1, because in this case Monty has no choice.

Applying Bayes's Theorem,
P(A|E) = p / (1+p)

In the canonical scenario, p=1/2, so P(A|E) = 1/3, which is the canonical solution. If p=0, P(A|E) = 0, so you can switch and win every time (when Monty opens B, that it). If p=1, P(A|E) = 1/2, so in that case it doesn't matter whether you stick or switch.

When Monty opens C, P(A|E) = (1-p) / (1+p).

7) And finally, here is a new problem I just came up with:

If you meet a man with (naturally) red hair, what is the probability that neither of his parents has red hair?

Hints: About 2% of the world population has red hair. You can assume that the alleles for red hair are purely recessive. Also, you can assume that the Red Hair Extinction theory is false, so you can apply the Hardy–Weinberg principle.

Solution to this one next week!

Please let me know if you have suggestions for more problems. An ideal problem should meet at least some of these criteria:

1) It should be based on a context that is realistic or at least interesting, and not too contrived.

2) It should make good use of Bayes's Theorem -- that is, it should be easier to solve with BT than without.

3) It should involve some real data, which the solver might have to find.

4) It might involve a trick, but should not be artificially hard.

If you send me something that is not under copyright, or is usable under fair use, I will include it in the next edition of Think Stats and add you to the contributors list.

My favorite Bayes's Theorem problems

2011-10-20T07:17:00.000-07:00

This week: some of my favorite problems involving Bayes's Theorem. Next week: solutions.

1) The first one is a warm-up problem. I got it from Wikipedia (but it's no longer there):

Suppose there are two full bowls of cookies. Bowl #1 has 10 chocolate chip and 30 plain cookies, while bowl #2 has 20 of each. Our friend Fred picks a bowl at random, and then picks a cookie at random. We may assume there is no reason to believe Fred treats one bowl differently from another, likewise for the cookies. The cookie turns out to be a plain one. How probable is it that Fred picked it out of Bowl #1?

This is a thinly disguised urn problem. It is simple enough to solve without Bayes's Theorem, but good for practice.

2) This one is also an urn problem, but a little trickier.

The blue M&M was introduced in 1995. Before then, the color mix in a bag of plain M&Ms was (30% Brown, 20% Yellow, 20% Red, 10% Green, 10% Orange, 10% Tan). Afterward it was (20% Blue , 20% Green, 16% Orange, 14% Yellow, 13% Red, 13% Brown).

A friend of mine has two bags of M&Ms, and he tells me that one is from 1994 and one from 1996. He won't tell me which is which, but he gives me one M&M from each bag. One is yellow and one is green. What is the probability that the yellow M&M came from the 1994 bag?

3) This one is from one of my favorite books, David MacKay's "Information Theory, Inference, and Learning Algorithms":

Elvis Presley had a twin brother who died at birth. What is the probability that Elvis was an identical twin?

Two people have left traces of their own blood at the scene of a crime. A suspect, Oliver, is tested and found to have type O blood. The blood groups of the two traces are found to be of type O (a common type in the local population, having frequency 60%) and of type AB (a rare type, with frequency 1%). Do these data (the blood types found at the scene) give evidence in favour [sic] of the proposition that Oliver was one of the two people whose blood was found at the scene?

5) I like this problem because it doesn't provide all of the information. You have to figure out what information is needed and go find it.

According to the CDC, ``Compared to nonsmokers, men who smoke are about 23 times more likely to develop lung cancer and women who smoke are about 13 times more likely.''
If you learn that a woman has been diagnosed with lung cancer, and you know nothing else about her, what is the probability that she is a smoker?

6) And finally, a mandatory Monty Hall Problem. First, here's the general description of the scenario, from Wikipedia:

Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say Door A [but the door is not opened], and the host, who knows what's behind the doors, opens another door, say Door B, which has a goat. He then says to you, "Do you want to pick Door C?" Is it to your advantage to switch your choice?

The answer depends on the behavior of the host if the car is behind Door A. In this case the host can open either B or C. Suppose he chooses B with probability p and C otherwise. What is the probability that the car is behind Door A (as a function of p)?

If you like this problem, you might also like the Blinky Monty Problem.

Solutions next week!

Please let me know if you have suggestions for more problems. An ideal problem should meet at least some of these criteria:

1) It should be based on a context that is realistic or at least interesting, not too contrived.

2) It should make good use of Bayes's Theorem -- that is, it should be easier to solve with BT than without.

3) It should involve some real data, which the solver might have to find.

4) It might involve a trick, but should not be artificially hard.

If you send me something that is not under copyright, or is usable under fair use, I will include it in the next edition of Think Stats and add you to the contributors list.

The Blinky Monty Problem

2011-10-17T13:02:00.000-07:00

I read Jason Rosenhouse's book about The Monty Hall Problem recently, and I use the problem as an example in my statistics class. Last semester I wrote a variation of the problem that turns out to be challenging, and a motivating problem for Bayesian estimation. Here's what I call the "Blinky Monty Problem."

Suppose you are on Let's Make a Deal and you are playing the Monty Hall Game, with one twist. Before you went on the show you analyzed tapes of previous shows and discovered that Monty has a tell: when the contestant picks the correct door, Monty is more likely to blink.

Of the 18 shows you watched, the contestant chose the correct door 5 times, and Monty blinked three of those times. Of the other 13 times, Monty blinked three times.

Assume that you choose Door A. Monty opens door B and blinks. What should you do, and what is your chance of winning?

To get started, let's review the standard version of the Monty Hall problem (from Wikipedia):

"Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say A [but the door is not opened], and the host, who knows what's behind the doors, opens another door, say B, which has a goat. He then says to you, "Do you want to pick door C?" Is it to your advantage to switch your choice?"

We can solve this problem using Bayes's Theorem. There are two relevant hypotheses:

A: the car is behind door A
C: the car is behind door C

Before Monty opens door B, the prior probability for both hypotheses is 1/3. Now the evidence is

E: Monty opened door B.

To compute the posterior probabilities, we need the likelihoods P(E|A) and P(E|C).

P(E|A) is the probability of the evidence if the car is behind door A. In this case Monty has a choice; he could have opened door B or C. In the canonical version of the problem, he chooses at random, so

P(E|A) = 1/2

P(E|C) is the probability of the evidence if the car is behind door C. In this case Monty has no choice, so

P(E|A) = 1

Applying Bayes's Theorem yields:

P(A|E) = 1/3
P(C|E) = 2/3

So if you switch to Door C, you increase the chances of winning from 1/3 to 2/3. This is the standard answer to the standard version of the problem. If you are not familiar with the Monty Hall problem and this answer comes as a shock to you, please don't take it out on me. Go read about it, starting here, and come back when you have calmed down.

Are you ready to go on? Ok. We can do a similar analysis for the Blinky Monty problem, but now we have more evidence to consider:

E: Monty opens Door B and blinks.

According to the scouting report, the probability that Monty blinks if the car is behind A is 3/5, so

P(E|A) = (1/2)(3/5) = 3/10

If the car is behind C, Monty blinks 3/13 of the time, so

P(E|C) = (1)(3/13) = 3/13

Plugging in Bayes' Theorem, we get P(A|E) = 0.51, so in this scenario you are slightly better off sticking with Door A.

But is that all there is to it? No! And here's where this blog earns the name "Probably Overthinking It." Remember that the probabilities for Monty's blinking are estimates based on a sample, so we have some uncertainty about them.

To take that uncertainty into account, we have to

1) Use the scouting report to estimate the distribution of P(blink|A) and P(blink|C).

2) For each pair of values from those distributions, compute P(A|E).

3) Apply the law of total probability to get the overall P(A|E).

That might sound complicated, but it is straightforward to estimate computationally.

First, let me take a minute to make fun of frequentists. According to conventional hypothesis testing, the data from the scouting report is not statistically significant. If the null hypothesis is that Monty has the same chance of blinking regardless of the location of the car, the p-value of the sample we saw is 0.27 (I used Fisher's exact test, computed using this online calculator). We can't reject the null hypothesis, so if we play by the rules of conventional hypothesis testing, I guess that means we can't take advantage of Monty's tell. If you are a committed frequentist, you should stop reading now.

Good. Now that the riff-raff are gone, let's proceed. Here's how we make the posterior distribution for P(blink|doorA):

doorA = MakeUniformSuite(0.0, 1.0, 101, name='Door A')
evidence = 3, 2
Update(doorA, evidence)

MakeUniformSuite() makes a Pmf of 101 values equally spaced between 0 and 1, with equal probability. Update() does the usual Bayesian update:

def Update(suite, evidence):

"""Updates a suite of hypotheses based on new evidence.

Modifies the suite directly; if you want to keep the original, make

a copy.

Args:

suite: Pmf object

evidence: whatever kind of object Likelihood expects

"""

for hypo in suite.Values():

likelihood = Likelihood(evidence, hypo)

suite.Mult(hypo, likelihood)

suite.Normalize()

def Likelihood(evidence, hypo):

"""Computes the likelihood of the evidence assuming the hypothesis is true.

Args:

evidence: a tuple of (number of heads, number of tails)

hypo: float probability of heads

Returns:

probability of tossing the given number of heads and tails with a

coin that has p probability of heads

"""

heads, tails = evidence

p = hypo

return pow(p, heads) * pow(1-p, tails)

If you are not familiar with that, you can read Chapter 8 of Think Stats, or see this blog post.

To find P(blink|C), it's pretty much the same:

doorC = MakeUniformSuite(0.0, 1.0, 101, name='Door C')

evidence = 3, 10

Update(doorC, evidence)

And to finish it off, we apply the law of total probability:

print TotalProbability(doorA, doorC, ProbWinning)

TotalProbability() enumerates the values from doorA and doorC, and calls ProbWinning() for each pair:

def TotalProbability(pmf1, pmf2, func):

"""Enumerates pairs from the Pmfs, calls the func, and returns

the total probability.

func: takes a value from each Pmf and returns a probability.

"""

total = 0.0

for x, px in pmf1.Items():

for y, py in pmf2.Items():

if px and py:

total += px * py * func(x, y)

return total

Finally, ProbWinning() takes P(blink|A) and P(blink|C) and returns P(A|E):

def ProbWinning(pbA, pbC):

"""Computes the probability that the car is behind door A:

pbA: probability that Monty blinks if the car is behind A

pbC: probability that Monty blinks if the car is behind C

"""

pea = 0.5 * pbA

pec = 1.0 * pbC

pae = pea / (pea + pec)

return pae

You can download all this code here. If you run it, you'll see the posterior distributions for P(blink|A) and P(blink|C):

The most likely values are 0.6 and 0.23, but since we don't have much data, the uncertainty is still large.

Taking it all together, the probability that the car is behind Door A, P(A|E), is 0.52. So you are slightly better off sticking with Door A.

To summarize:

1) In the canonical version of the game, P(A|E) = 1/3, so you are better off switching.

2) In the Blinky Monty version, if we take P(blink|A) and P(blink|C) as givens, we estimate P(A|E) = 0.51, so you are slightly better off sticking with Door A.

3) If you are a frequentist, you conclude that Monty's tell is not statistically significant, so you ignore it, switch, and lose most of the time.

4) If you are a Bayesian, you can take advantage of the tell and maximize your chance of winning. If Monty blinks and you stick, you win 52% of the time. If he doesn't blink and you switch, you win 87% of the time. Since he blinks 40% of the time, overall you can expect to win 73% of the time. So the tell provides a substantial advantage.

Exercise for the reader: download my code and modify it to analyze the case where Monty doesn't blink. Hint: you only have to change two lines.

EDIT October 19, 2011: I got a comment on reddit that raises a legitimate concern about data collection. In the scenario I presented, if you review tapes of the show and look for several potential tells (blinking, facial tics, verbal habits, etc.) and you choose one that seems to be correlated with the outcome, you could easily get fooled (see this post on repeated tests). This is a fair warning for both the frequentist and Bayesian analyses.

So to be more careful, maybe I should say that you were tipped off to the possibility that Monty's tell is a blink before you reviewed the tape, and that's the only data you collected.

EDIT January 26, 2012: In the comments below, Jacob pointed out an error in my calculation. As a result, I went back and changed the scenario to match my calculations, which forced me to make a few other revisions. I think it is all settled now, but let me know if you disagree.

Repeated tests: how bad can it be?

2011-10-06T05:03:00.000-07:00

A couple of weeks ago Red Dave wrote this blog post about the dangers of A/B testing with repeated tests. I think he's right, and does a good job of explaining why. But to get a better understanding of the issue, I wrote some simulations, and things can be even worse than Dave suggests.

As an example of A/B testing, suppose you are advertising a business and you develop two versions of an ad, A and B, with different text, graphics, etc. You want to know which ad is more effective.

If you post your ads online, you can find out how many times each version was shown (number of impressions), and how many people clicked on each ad (number of clicks). The number of clicks per impression is the "click-through rate," or CTR. At least, that's what Google calls it.

If you show each ad 1000 times, you might see a CTR of 1% for Version A and 2% for Version B. To see whether this difference is significant, you can use a chi-square test, which compares the actual number of clicks for each ad with the expected number. The result is a p-value, which tells you the chance of seeing a difference as big as that if the ads are actually the same. If the p-value is small, the effect is unlikely to be due to chance, so it is more likely to be real.

So far, this is just classical hypothesis testing. No problem. Here's where things get ugly.

If, instead of showing each ad 1000 times, you run the ad for a while, and then check for a statistically-significant difference, run a while longer, check again, and repeat until you get a significant difference, then what you are doing is Bogus.

And when I say "Bogus", I mean capital-B Bogus, because the p-values you compute are not going to be off by a factor of 2, which would actually be acceptable for most purposes. They are not going to be off by a factor of 10, which would make them meaningless. They are going to be completely and totally wrong, to the point where, if you run long enough, you are nearly certain to see a "statistically significant" difference, even if both versions are the same.

To see why, let's look at the case where you check for significance after every ad impression. Suppose you show version A or version B with equal probability, and in either case the click-through rate is 50%. After 100 impressions, you expect to show A and B 50 times each, and you expect 25 clicks on each.

After each impression, we can compute the CTR for A and B, and plot the difference as a random walk. Here's what it looks like for 10 simulations of 200 impressions:

When the number of trials is small, the difference in CTR can be big, but for n>50 the law of large numbers kicks in and the difference gets small.

Instead of computing the difference, we can compute the chi-square statistic, which yields another random walk. Here's what it looks like after 1000 steps:

After an initial transient, the distribution of the lines stabilizes, which shows that for large N, the distribution of the chi-square stat does not depend on N, which is exactly why we can use the chi-square distribution to compute p-values analytically.

SciPy provides a function that evaluates the chi-squared CDF. The following function uses it to look up chi-squared values and return the corresponding p-value:

import scipy.stats

def Pvalue(chi2, df):
"""Probability of getting chi2 or greater from a chi-squared distribution.

chi2: observed chi-squared statistic
df: degrees of freedom
"""
return 1 - scipy.stats.chi2.cdf(chi2, df)

In this example, the number of degrees of freedom is 3, because there are four values in the table (hits and misses for A and B), but they are constrained to add up to N.

Now here's the key question: if we take values that are actually from a chi-squared distribution and look up the corresponding p-values, what is the distribution of the resulting p-values?

Answer: they are uniformly distributed. Why? Because we are using the chi-squared CDF to get their percentile ranks. By the definition of the CDF, 1% of them fall below the 1st percentile, 2% of them fall below the 2nd percentile, and so on.

If we plot the p-values as a random walk, here's what it looks like (20 simulations, 2000 trials):

For large N, the lines are spread out uniformly from 0 to 1, as expected.

The red line shows the significance criterion, alpha, which determines the false positive rate. If alpha=5%, then at any point in time, we expect 5% of the lines to be below this threshold. And in the figure, that looks about right.

If we stop the test at a pre-determined time, or at a randomly-chosen time, we expect to find one line out of twenty below the threshold. But if you do repeated tests, checking for significance after every trial, then the relevant question is how many of these lines cross below the threshold, ever?

In this example, it is almost 50%. After 10000 trials, it is 85%. As the number of trials increases, the false positive rate approaches 100%.

A standard quip among statisticians goes like this: "A frequentist is a person whose long-run ambition is to be wrong 5% of the time." Here's an update: "A/B testing with repeated tests is totally legitimate, if your long-run ambition is to be wrong 100% of the time."

You can download the code I used to generate these figures from thinkstats.com/khan.py. The file is named after the Khan Academy because that's the context of the blog post that brought the problem to my attention.

Note: in these examples, I used a CTR of 50%, which is much higher than you are likely to get with an online ad campaign: 1% is considered good. If you run these simulations with more realistic CTRs, the random walks get jaggier, due to the statistics of small numbers, but the overall effect is the same.

The Jimmy Nut Company problem

2011-08-29T18:33:00.000-07:00

I read with interest an article by Robert Hayden called “Advice to mathematics teachers on evaluating introductory statistics textbooks.” It makes two claims:

Statistics textbooks should emphasize real problems with real data. I couldn't agree more.
You should choose a statistics textbook written by a statistician at a large university who has good research publications. I couldn't agree less.

On the second point, Prof. Hayden argues that statistics books contain more errors than other math books, and that the errors are more often conceptual, not just typos. The reason, he explains, is that “introductory statistics textbooks are often written by people with little or no training in statistics.”

As an example (in fact, the only example) he cites this problem from an unnamed popular textbook:

“Jimmy Nut Company advertises that their nut mix contains 40% cashews, 15% brazil nuts, 20% almonds, and only 25% peanuts. The truth in advertising investigators took a random sample (of size 20 lb) of the nut mix and found the distribution to be as follows:

Cashews Brazil Nuts Almonds Peanuts

6 lb 3 lb 5 lb 6 lb

At the 0.01 level of significance, is the claim made by Jimmy Nuts true?”

Prof. Hayden hates this problem. He thinks it is ill-posed, and evidence of everything wrong with statistics books. He explains: “The problem here is that the chi-squared goodness of fit test applies only to categorical (discrete) data. It compares the actually observed counts in each category to the counts we would expect if the hypothesis being tested were true.”

This is true, and if you blindly plug the numbers from the problem into the formula for a chi-squared statistic, the result is meaningless. In fact, if you express the given data in other units, you can make the resulting p-value as big or small as you like.

But that doesn’t mean the textbook author has made a “gross error,” because the question doesn’t say “plug these numbers in blindly and write down the result without thinking.” The question doesn’t even specify the chi-squared statistic.

The question, as written, is perfectly sensible, moderately challenging, and interesting. It requires thought and care to formulate the question in a way that leads to a meaningful answer, and it is exactly the kind of question someone doing statistics is likely to be asked.

I don’t know what the textbook author had in mind, but there is nothing wrong with the question as written. And to prove it, I’m going to answer it. I follow the procedure I recommend in my article "There is only one test."

The first step is to make a model of the null hypothesis, which is that the proportion of nuts is as advertised. But if the proportions are always as advertised, the data observed would be impossible. So we need a random model that characterizes the variation we expect to see in a sample.

What are the sources of variation in a 20 lb sample of nuts? One source is mechanical error; that is, the people and machines that weigh and mix nuts might be miscalibrated. If the accuracy of this process varies from day to day, we could treat it as a random variation. To quantify this effect, we would need some information about how the nuts are processed.

Another source of error, and probably the one the textbook author had in mind, is discrete sampling. If we take a small sample, we would not expect the proportions to be exactly as advertised. To model this kind of variation, we have to count individual nuts.

To convert the data from weights to counts, we have to know how much each kind of nut weighs. This step requires some research, which I think is an interesting element of the problem. Since we live in a post-Google information utopia, it is reasonable to expect students to answer this question, and useful for them to practice.

I found several sources with information about the weight of typical nuts, some more authoritative than others. The one I use is “Nuts for Nutrition,” by Alice Henneman, MS, RD, Extension Educator at the University of Nebraska Cooperative Extension - Lancaster County. In addition to the credentials and affiliation of the author, and the witty article title, I chose this source because it provides numbers for all four nuts in one table, so even if the methodology that produced them is not ideal, it is probably consistent, which means that the relative values are likely to be good even if the absolute values are not. From the table:

Nuts per ounce: Cashew 16-18, Brazil nut 6-8, Almonds 20-24, Peanuts 28.

To start I use the middle of each range. With this data we can convert observed values from pounds to counts and convert the expected proportions from “by weight” to “by count”. Here’s a function that does it:

def ConvertToCount(sample, count_per):
   """Convert from weight to count.

   sample: Hist that maps from category to weight in pounds
   count_per: dict that maps from category to count per ounce
   """
   for value, count in sample.Items():
       sample.Mult(value, 16 * count_per[value])

The code examples here use the Pmf and Cdf libraries from Think Stats. Here’s the code that processes the information we’re given:

   sample = dict(cashew=6, brazil=3, almond=5, peanut=6)
   count_per = dict(cashew=17, brazil=7,
almond=22, peanut=28)

   observed = Pmf.MakeHistFromDict(sample)
   ConvertToCount(observed, count_per)

   advertised = dict(cashew=40, brazil=15,
almond=20, peanut=25)
   expected = Pmf.MakePmfFromDict(advertised)
   ConvertToCount(expected, count_per)
   expected.Normalize(observed.Total())

Here are the results in tabular form:

Nut       Expected     Observed    % error    price per pound
cashew    2266         1632        -28 %      $10
brazil     350          336        - 4 %      $8
almond    1467         1760        +20 %      $8
peanut    2333         2688        +15 %     $3

So there are 28% fewer cashews than expected, and 20% more almonds. Is this suspicious? Intuitively, this looks like more deviation than I expect in a sample this large. And the relationship between “error” and price makes it look even worse.

Let’s make that more rigorous. The next step to choose a test statistic. Mean relative error would be a reasonable choice. Or I could devise a statistic that measures the excess of cheap nuts and lack of expensive nuts. But to keep it simple I'll use the chi-square statistic. The primary advantage of chi-square is that we can compute p-values analytically, which is efficient. But since I am planning to use simulation, the only advantage is that it is conventional.

For the given data the chi-square statistic is 291, which doesn’t mean anything. To make it mean something, we have to compute the p-value, which is the probability of seeing a chi-square statistic so large under the null hypothesis.

Which means we need to get back to our model of the null hypothesis. One option is to imagine choosing nuts one at a time from a vat that contains a large number of nuts mixed in the advertised proportions. Here’s what that looks like:

   num_nuts = observed.Total()
   cdf = Cdf.MakeCdfFromHist(expected)
   t = cdf.Sample(num_nuts)
   hist = Pmf.MakeHistFromList(t)
   chi2 = ChiSquared(expected, simulated)

If we run that 1000 times, we can count how many samples yield chi2 > 290. And the answer is...0. In fact, after 1000 simulations, the largest test statistic is 20. So the p-value is much smaller than 1/1000 and it is unlikely that the observed sample came from the advertised mixture, at least under this model.

Since there was some uncertainty in the data we collected, we should do a perturbation analysis. For example, Jimmy might claim that he used big nuts, which would make the counts smaller and the results less significant. So let’s give Jimmy the benefit of the doubt and run the numbers with these counts per ounce:

Nuts per ounce: Cashew 12, Brazil nut 5, Almonds 18, Peanuts 24.

With these data, the observed test statistic is only 215 (down from 290). But the chance of seeing a value that large under the null hypothesis is still negligible.

More credibly, Jimmy might claim that the model of the null hypothesis is unrealistic because the nuts are not perfectly mixed, so the sampling process is not independent. In fact, correlation between successive draws would increase variation in the sample. To see why, imagine the extreme scenario where the nuts are not mixed at all; in that case, most samples contain only one kind of nut.

To model this scenario, we can make a list that represents the vat of nuts the sample is drawn from. Initially the vat is completely unmixed. Then we "stir" by choosing random pairs and swapping them.

def MakeVat(expected, num_nuts, factor=10, stir=0.0):
"""Makes a list of nuts with partial stirring."""
t = []
for value, freq in expected.Items():
t.extend([value] * int(freq * factor))

[RandomSwap(t) for i in xrange(int(num_nuts*factor*stir))]

return t

expected is a Pmf that represents the expected distribution. num_nuts is the total number of nuts in the sample. factor is the size of the vat, in multiples of num_nuts. And stir is the number of swaps to perform per nut in the vat.

I set factor to 10, which means nuts are mixed in 200 lb batches. With stir = 0, the nuts are unmixed and the p-value is 1; that is, it is certain that we will see a large deviation from the advertised proportions. If stir > 2, the nuts are well mixed and the p-value is 0, as we saw before.

Between 1.1 and 1.3 swaps per nut, the p-value drops quickly from 1 to 0. So Jimmy's argument is plausible: if the nuts are not well mixed, that could explain the deviation we saw in the sample.

In summary, what can we conclude about the Jimmy Nut Company?

1) If the sample is drawn from a well mixed batch of nuts, it provides strong evidence that the mix is not as advertised.

2) But if the nuts are not well mixed, that might explain the deviation we saw. We could test that hypothesis by mixing the vat more thoroughly or taking small samples from several parts of the vat.

And what can we conclude about the Jimmy Nut Company problem?

1) I think it's an excellent problem, as written. It requires careful thinking and a little research, and the results could lead to interesting class discussion.

2) I see no evidence that this problem indicates a conceptual error on the part of the textbook writer.

In summary, I agree with Prof. Hayden that statistics students should work with real data, but respectfully disagree with his assertion that statistics books should be written by statisticians. I think they should be written by computer scientists.

The code I used for this article is here.

UPDATE August 30, 2011: With the help of Google Books, I found the source of this problem: Understandable statistics: concepts and methods by Charles Henry Brase and Corrinne Pellillo Brase. The problem appears in the 5th edition, but not the current 9th edition.

Upcoming webcast: Only One Test

2011-08-16T07:27:00.000-07:00

The nice people at O'Reilly Media have asked me to do a webcast or two on topics related to Think Stats. The first one is scheduled for October 4 at 1pm EDT. Here is the announcement.

There's Only One Test

Tuesday, October 4, 2011

Presented by: Allen B. Downey

Duration: Approximately 60 minutes.

Cost: Free

Can I use a t-test if my data are non-normal? What if the sample size is small? What's an exact test? And how many tests are there, anyway?

People working with real data are often confused about hypothesis testing and paralyzed by the number of tests and their requirements. But for many common problems, you don't need specialized tests at all.

In this talk, I present a framework for using simple simulations to estimate p-values. This framework works well for a wide range of problems, which is why I say that there's only one test.

You don't need prior knowledge of statistics (and you might be better off without it). I show programming examples in Python, but the talk is accessible to anyone with basic programming skills.

About Allen B. Downey

Allen B. Downey, author of Think Stats. Allen is a Professor of Computer Science at Olin College of Engineering. He has taught at Wellesley College, Colby College and U.C. Berkeley. He has a Ph.D. in Computer Science from U.C. Berkeley and Master's and Bachelor's degrees from MIT.

Questions? Please send email to webcast@oreilly.com

Here are the slides I used for this webcast:

EDIT Oct 12, 2011: Posted final version of slides.

A hierarchical Bayesian model of pond scum

2011-06-08T07:23:00.000-07:00

This week I am working with one of my colleagues, the extraordinary biologist Jean Huang, on an interesting problem related to bioinformatics. This project is my first attempt at implementing a hierarchical Bayesian model in Python, and it has been a pleasant surprise. If you are familiar with tree-processing algorithms, it's pretty straightforward. First the likelihoods propagate up the tree, then the updates propagate down.

Here is the background: Prof Huang and her students collect samples of pond water and culture the bacteria they find under narrow-spectrum artificial light, trying to find species that photosynthesize at different wavelengths.

To identity the species in each culture, they amplify the 16S ribosomal RNA gene, send it out to be sequenced, and then look up the sequence in a database like greengenes. The result is either the name of a species that has been cultured, the genus of an uncultured species, or an unknown genus.

They usually start by identifying 15 samples from each culture. For example, one of their cultures yielded 9 samples of one species, three of another, and one each of three more.

Looking at the results for a given culture, biologists would like to be able to answer these questions:

1) How many more species are there likely to be?
2) What is the likely prevalence of each species?
3) If we test m additional samples, how many new species are we likely to find?
4) Given a limited budget, which cultures warrant additional sampling?

To answer these questions, I implemented a hierarchical Bayesian model where the top level is "How many species are there?" and the second level is "Given that there are n species, what is the prevalence of each species?" The prevalence of a species is its proportion of the population.

Knowing how many species there are helps model their prevalences. In the simplest case, if we know that there is only 1 species, the prevalence of that species must be 100%. If we know there are 2 species, I assume that the distribution of prevalences for both species is uniform from 0 to 100, so the expected prevalence for both is 50%, which makes sense.

If we know there are 3 species, it is less obvious what the prior distribution should be; the Beta distribution is a natural choice because

1) If we know there are k species, we can use Beta(1, k-1) as a prior, and the expected value is 1/k, which makes sense. For k=2, the result is a uniform distribution, so that makes sense, too.

2) For this scenario, the Beta distribution is a conjugate prior, which means that after an update the posterior is also a Beta distribution.

We can represent a Beta distribution with an object that keeps track of the parameters:

class Beta:
def __init__(self, yes, no):
"""Initializes a Beta distribution.
yes and no are the number of
successful/unsuccessful trials.
"""
self.alpha = yes+1
self.beta = no+1

Update couldn't be easier.

def Update(self, yes, no):
"""Updates a Beta distribution."""
self.alpha += yes
self.beta += no

At the top level ("How many species are there?") the prior is a uniform distribution from 1 to 20, at least for now. I could replace that with a distribution that reflects more of the domain knowledge biologists have about the experiment; for example, based on how the cultures are processed, it would be rare to find more than 10 species with any substantial prevalence.

Now that the priors are in place, how do we update? Here's the top level:

def Update(self, evidence):
"""Updates based on observing a given taxon."""
for hypo, prob in self.hypos.Items():
likelihood = hypo.Likelihood(evidence)
if likelihood:
self.hypos.Mult(hypo, likelihood)
else:
# if a hypothesis has been ruled out, remove it
self.hypos.Remove(hypo)

self.hypos.Normalize()

The evidence is a string that indicates which taxon (species or genus) was observed. self.hypos is a Pmf that maps from a hypo to its probability (Pmf is provided by the thinkstats package, a collection of modules used in my book, Think Stats). Each hypo is a hypothesis about how many taxons there are. We ask each hypo to compute the likelihood of the evidence, then update hypos accordingly. If we see N different taxons, all hypotheses with k<N are eliminated.

Having updated the top level, we pass the evidence down to the lower level and update each hypothesis.

# update the hypotheses

for hypo, prob in self.hypos.Items():

hypo.Update(evidence)

To compute the likelihood of the evidence under a hypothesis, we ask "What is the probability of observing this taxon, given our current belief about the prevalence of each taxon?" There are two cases:

1) If the observed taxon is one we have seen before, we just look up the current Beta distribution and compute its mean:

def Mean(self):

"""Computes the mean of a Beta distribution."""

return self.alpha / (self.alpha + self.beta)

2) If we are seeing a taxon for the first time, we get the likelihoods for all unseen taxons and add them up.

To update the lower level distributions, we loop through the Beta distributions that represent the prevalences, and update them with either a hit or a miss.

def Update(self, evidence):

"""Updates each Beta dist with new evidence.

evidence is a taxon

"""

for taxon, dist in self.taxa.iteritems():

if taxon == evidence:

dist.Update(1, 0)

else:

dist.Update(0, 1)

And that's it. After doing a series of updates, the top level is the posterior distribution of the number of taxons, and each hypothesis at the lower level is a set of Beta distributions that represent the prevalences.

Let's look at an example. Suppose we see 11 of one taxon, 4 of another, and 2 of a third. The posterior distribution of the number of taxons is:

The probability is 47% that there are only 3 taxons, 28% that there are 4, and 13% that there are 5.

By summing over the hypotheses, we can generate distributions for the prevalence of each taxon:

The median prevalence of taxon A is 58%; the credible interval is (39%, 75%). For taxons B and C the medians are 23% and 13%. For taxons D, E, and F, there is a high probability that the prevalence is 0, because these taxons have not been observed.

To predict the number of new taxons we might find by sequencing additional samples, I use Monte Carlo simulation. Here is the kernel of the algorithm:

for i in range(m):

taxon = self.GenerateTaxon()

taxons.append(taxon)

meta.Update(taxon)

curve = MakeCurve(taxons)

m is the number of samples to simulate. Each time through the loop, we generate a random taxon and then update the hierarchy. The result is a rarefaction curve, which is the number of taxons we observe as a function of the number of samples.

To generate a random taxon, we take advantage of the recursive structure of the hierarchy; that is, we ask each hypothesis to choose a random taxon, and then choose among them in accordance with the probability for each hypothesis:

def GenerateTaxon(self):

"""Chooses a random taxon."""

pmf = Pmf.Pmf()

for hypo, prob in sorted(self.hypos.Items()):

taxon = hypo.GenerateTaxon()

pmf.Incr(taxon, prob)

return pmf.Random()

The algorithm the hypotheses use to generate taxons is ugly, so I will spare you the details, but ultimately it depends on the ability to generate Beta variates, which is provided in Python's random module.

Given a sequence of taxons (real or simulated), it is easy to generate a rarefaction curve:

def MakeCurve(sample):

"""Makes a rarefaction curve for the given sample."""

s = set()

curve = []

for i, taxon in enumerate(sample):

s.add(taxon)

curve.append((i+1, len(s)))

return curve

To show what additional samples might yield, I generate 100 rarefaction curves and plot them:

The first 17 points are random permutations of the observed data, so after 17 samples, we have seen 3 taxa on every curve. The last 15 points are based on simulated data. The lines are shifted slightly so they don't overlap; that way you can eyeball high-probability paths.

Even after m=15 additional samples, there is a good chance that we don't see any more taxons -- but there is a reasonable chance of seeing 1-2, and a small chance of seeing 3.

Finally, we estimate the probability of seeing additional taxons as a function of the number of addition samples:

With 5 additional samples, the chance of seeing at least one additional taxon is 24%; with 10 samples it's 40%, and with 15 it's 53%.

These estimates should help biologists allocate their budget for additional samples is order to minimize the chance of missing an important low-prevalence taxon.

Next steps: I am working with Prof. Huang and her students to report the results most usefully; then we are thinking about distributing the code or deploying it as a web app.

More hypotheses, less trivia

2011-06-02T15:18:00.000-07:00

In my last post I wrote about hypothesis testing with simulation, and one of the readers at reddit/r/statistics wrote this comment:

Is there a book or similar reference [with] examples of building such simulations? I never have enough confidence to run such simulations and I guess that is mainly because I have only seen them for trivial problems. A book or a tutorial with a number of problems with various complexity would help a lot.

My answer to the first part (is there a good book?) is "I don't know." My answer to the second part (non-trivial examples) is this post with some less trivial examples.

Difference in means

If you think there is a difference between two groups, the most common test is for a difference in the means. In Think Stats, I apply this test to pregnancy length and birth weight for first babies and others.

In this case the test statistic is obvious: compute the mean for each group and subtract. The only hard part is deciding whether to do a one-sided test or a two-sided test. It depends on what hypothesis you are testing. If you think men are taller than women, or that the new drug is better than the old drug, then a one-sided test is appropriate. If you think there is a difference between the groups but you don't know what it is, you can use a two-sided test; that is, use the absolute value of the difference as the test statistic.

This decision can seem arbitrary, but don't worry. The effect on the p-value is just a factor of two, and (as I argued last time) we only care about the order of magnitude. Whether the p-value is 2% or 4% or 8% doesn't really matter.

The null hypothesis is that the two groups are actually the same. There are several ways to implement this model:

Parametric: Merge the groups. If the pooled distribution is well modeled by an analytic distribution, estimate the parameters (for example, the mean and variance of a normal distribution) and then generate random samples from the fitted distribution.

Resampling: Again, merge the groups, and then draw samples (with replacement) from the pooled empirical distribution.

Permutation: Or shuffle the sample and assign the elements to the groups at random. This is equivalent to resampling without replacement.

In this example, the difference in birthweight is significant (first babies are lighter) and the difference in pregnancy length is not. But it turns out that first babies are more likely to be early and late (and less likely to be on time). To test that hypothesis, I used the procedure I described last time for the casino problem.

Slope of a fitted line

In this post I fit a line to a time series and compute the slope. The test statistic is the estimated slope, so that part is easy.

The null hypothesis is that the apparent slope is due to chance, so I used a permutation model; that is, I shuffled the time series 1000 times and computed the slope each time.

Since my hypothesis is specifically that the slope is positive, I could have used a one-sided test, but to be conservative I reported a two sided p-value:

The red line shows the linear least squares fit, with slope 0.033; the p-value (chance of seeing an absolute slope as big as that) is 0.006, so you can either reject the null hypothesis or update your subjective degree of belief accordingly.

Notice that I didn't say anything about how I computed the p-value. Depending on where you publish, you might or might not have to be more specific.

Tail curvature

In this paper, I was characterizing the tail behavior of long-tailed distributions and trying to quantify how well Pareto and lognormal models captured that behavior.

If you plot the complementary CDF of a Pareto distribution on a log-log scale, you get a straight line (example here); for a lognormal distribution, the line drops off with increasing slope. To distinguish these behaviors, I wanted a robust measure of tail curvature, so I computed the first derivative numerically and the estimated the slope. Since slope of a derivative is the second derivative, I called the result "tail curvature" and used it as the test statistic.

In this example I had two hypotheses (HP for Pareto and HL for lognormal) and no null hypothesis. To compute P(E | HP) I chose a Pareto distribution that was a good fit for the tail of the observed distribution, then generated random samples, and computed tail curvature. Likewise for P(E | HL).

Here's how I reported the results:

For these cases, the curvature test provides some insight. The tail curvature of the Calgary dataset is 0.141, which has a negligible p-value under the Pareto model. This means that we can reject the hypothesis that the data are a sample from a Pareto distribution. The tail curvature of the ClarkNet dataset is 0.023, which has p > 0.5, which means we cannot reject the hypothesis that the sample comes from a long-tailed distribution.

I am not crazy about the "reject the null hypothesis" way of talking, but it was appropriate for the venue.

ANOVA

If you collect a sample of adult human heights in centimeters, the variance might be 94 cm^2. If you divide the sample into men and women, the variance in each group is smaller (49 and 44 cm^2), which suggests that there is a relationship between height and gender. To test whether that is significant, you can divide the sample at random and see how often the variance drops as much. That's an analysis of variance (ANOVA) test.

The nice thing about ANOVA is that is generalizes easily to more than two groups.

To summarize, the test statistic is the decrease in variance when you divide the sample into groups, also known as "assigning labels." The null hypothesis is that the mean is the same for all groups, so the simulation is to assign the labels at random.

Time of day effect?

Here's a recent post from reddit/r/statistics:

A reviewer wants me to test if various behavioural assays I performed were influenced by the time of day I performed them (i.e. biased).The assays have either binary or continuous responses. What is the best approach to this? It is as simple as fitting a GLM, response~time, and testing for a non-zero slope? That seems overly simplistic to me. I've also tried GAMs with gaussian or binomial error functions, but it doesn't seem to be doing what I expect it to do...Any thoughts would be greatly appreciated. I'm sure its stupidly simple, but I'm having trouble wrapping my head around the problem!

At first I suggested a permutation test to check whether time of day had a substantial effect. But the submitter added these details:

We assayed the behaviour of six species of spiders. The purpose of the study was to compare the behaviour of some species to others. For various reasons, we could not balance the time of day at which we assayed the behaviours for each species (the wonders of field biology...). It may be the case that these spiders have some sort of daily behavioural cycle, e.g. more active at night (however it could be the opposite, or really any function of time; we have no idea). I think the reviewer's concern is that this could influence the comparisons, e.g. if we assayed species 1 more often at night than species 2?
I suppose the explicit hypothesis I am interested in testing is that time of day has a significant effect on behavioural response, but I am unsure how to construct a model to test this, given that I have no idea if the relationship will be linear or some complex function.

So there are several questions on the table here:

1) Are the observed differences between the species significant?

I assume the answer is yes, or the submitter would not be trying to publish (see publication bias).

2) Is there a time of day effect?

One option is to divide the day into intervals, maybe 6 periods of 4 hours. Then for the continuous response variable, you could use an ANOVA test. For the binary variable you could use a logistic regression with one dummy variable for each interval.

Or if the question is whether there are any times of day that differ from the others, then an appropriate test statistic would be the maximum difference between any interval and the overall average (the next case study talks about this option in more detail).

3) If there is a time of day effect, can we control for it?

If it turns out that there is a time of day effect, it may or may not be possible to control for it. For example, if you only observed Species 1 in the morning and Species 2 in the afternoon, then there is no way to separate the effect of the two variables. In that case, the problem is experimental design, not statistics.

As a quick check, you could make a table with one row per species and one column per time interval, and count how many observations fall in each cell. If there are a lot of empty (or nearly empty) cells, you are in trouble.

If you make the intervals bigger, you will have more observations per cell, but you might obscure an actual time of day effect.

4) If we control for time of day, are the differences between the species still significant?

If you get this far, you could do ANOVA or regression with dummy variables for each species and for each interval. This is what some people suggested on reddit (if I understood correctly), but I would suggest spending some time on steps (2) and (3) before diving into (4).

Testing pair-wise differences

Here's another recent post from reddit/r/statistics:

My data can be divided into 4 samples. I wanted to test whether the distribution of a binary variable differs significantly between the samples, so I did 6 t-tests - one for each pair of samples.
However, I also need to analyze how another categorical variable differs between the samples, this time with 3 possible values.
Is it okay to compute new t-tests, each time for a pair of values and a pair of samples, for a total of 6*3 tests? The problem is that this is too detailed: I don't want to test each pair of values, only each pair of samples.
SAS also gives me the option (in proc freq) of calculating a chi-square test for all 4 samples, but this is not enough information: it only tells me that not all of the means are equal.
And then there's a one-way ANOVA test, which can be used to compare the means of >2 samples, but apparently can't be used for categorical variables with >2 values. Or can it?
How would you analyze this data?

If you run lots of tests, some of them are likely to be "significant" even if there is no real effect going on. The problem is demonstrated in this xkcd cartoon. So in general I think it is a good idea to apply one big test to the entire dataset, not lots of little tests.

In this case, the question seems to be whether there are significant differences between any pair of groups. So an appropriate test statistic would be the maximum difference between any pair.

The null hypothesis is that the groups are the same, so an easy implementation is to divide the sample into groups at random.

For categorical variables, you can think of the distribution within each group as a vector and use a distance function to compute the difference between groups. If the groups are roughly the same size, you could use a 1-norm or a 2-norm. If not, you should probably normalize the vectors to length 1 before computing norms.

If you are not really interested in pair-wise differences, but you want to know whether any group deviates significantly from the expected, you could pool all the groups to get the expected distribution of the categorical variable, then compute a chi-squared statistic for each group (how much it deviates from the expected). Again the test statistic would be the maximum of the chi-square statistic.

The important thing here is that the choice of test statistic depends on the question you are asking. By using the maximum difference, you are effectively asking, "What is the probability of finding any group with this much deviance, just by chance?"

Running a series of tests

Here's one more question from reddit/r/statistics

I'm running a series of chi-square analyses and for several of them, I'm getting a message which states that "# cells have an expected count less than 5". If the results are significant, are they still worth reporting in spite of this message? I'm comparing participant selections of one of three alternatives when assigned to one of four conditions, examining pairs of comparisons.

The problem here is that the analytic version of the chi-squared test uses an approximation that is only accurate if the sample size is large enough; otherwise you have to use a different test.

But this is a perfect example of what I talked about in my last post: if you spend all your time in statistics class learning one damn test after another, and your entire cognitive capacity is consumed trying to remember when each of them is applicable, you have no time and no available brain cells to do the real work, which is thinking carefully about the what you are doing. [Note: I mean this as a general warning, not a criticism of the submitter.]

In this case the biggest hazard is not the sample size; it's "running a series of chi-square analyses." When I hear "series of tests," alarm bells go off in my head. If you don't know why, please review this xkcd cartoon.

So for this question, I have two suggestions:

1) If you use simulation, you don't have to worry about the expected count in each cell. If the sample size is small, that tends to make the p-value bigger, but that happens without any additional effort.

2) If you are running a series of tests, you should either (a) adjust the significance criterion using something like the Holm–Bonferroni method, or (b) treat the whole series as a single test, using the maximum of the chi-square values as a test statistic. I think (b) is simpler, and does a better job of capturing what you are trying to test.

That's all for now!

There is only one test!

2011-05-31T18:28:00.000-07:00

A few weeks ago I started reading reddit/r/statistics regularly. I post links to this blog there and the readers have given me some good feedback. And I have answered a few questions.

One of the most commons question I see is something like, "I have some data. Which test should I use?" When I see this question, I always have two thoughts:

1) I hate the way statistics is taught, because it gives students the impression that hypothesis testing is all about choosing the right test, and if you get it wrong the statistics wonks will yell at you, and

2) There is only one test!

Let me explain. All tests try to answer the same question: "Is the apparent effect real, or is it due to chance?" To answer that question, we formulate two hypotheses: the null hypothesis, H0, is a model of the system if the effect is due to chance; the alternate hypothesis, HA, is a model where the effect is real.

Ideally we should compute the probability of seeing the effect (E) under both hypotheses; that is P(E | H0) and P(E | HA). But formulating HA is not always easy, so in conventional hypothesis testing, we just compute P(E | H0), which is the p-value. If the p-value is small, we conclude that the effect is unlikely to have occurred by chance, which suggests that it is real.

That's hypothesis testing. All of the so-called tests you learn in statistics class are just ways to compute p-values efficiently. When computation was expensive, these shortcuts were important, but now that computation is virtually free, they are not.

And the shortcuts are often based on simplifying assumptions and approximations. If you violate the assumptions, the results can be misleading, which is why statistics classes are filled with dire warnings and students end up paralyzed with fear.

Fortunately, there is a simple alternative: simulation. If you can construct a model of the null hypothesis, you can estimate p-values just by counting. This figure shows the structure of the simulation.

The first step is to get data from the system and compute a test statistic. The result is some measure of the size of the effect, which I call "delta". For example, if you are comparing the mean of two groups, delta is the difference in the means. If you are comparing actual values with expected values, delta is a chi-squared statistic (described below) or some other measure of the distance between the observed and expected values.

The null hypothesis is a model of the system under the assumption that the apparent effect is due to chance. For example, if you observed a difference in means between two groups, the null hypothesis is that the two groups are the same.

The next step is to use the model to generate simulated data that has the sample size as the actual data. Then apply the test statistic to the simulated data.

The last step is the easiest; all you have to do is count how many times the test statistic for the simulated data exceeds delta. That's the p-value.

As an example, let's do the Casino problem from Think Stats:

Suppose you run a casino and you suspect that a customer has replaced a die provided by the casino with a ``crooked die''; that is, one that has been tampered with to make one of the faces more likely to come up than the others. You apprehend the alleged cheater and confiscate the die, but now you have to prove that it is crooked. You roll the die 60 times and get the following results:

Value 1 2 3 4 5 6

Frequency 8 9 19 6 8 10

What is the probability of seeing results like this by chance?

The null hypothesis is that the die is fair. Under that assumption, the expected frequency for each value is 10, so the frequencies 8, 9 and 10 are not surprising, but 6 is a little funny and 19 is suspicious.

To compute a p-value, we have to choose a test statistic that measures how unexpected these results are. The chi-squared statistic is a reasonable choice: for each value we compare the expected frequency, exp, and the observed frequency, obs, and compute the sum of the squared relative differences:

def ChiSquared(expected, observed):
total = 0.0
for x, exp in expected.Items():
obs = observed.Freq(x)
total += (obs - exp)**2 / exp
return total

Why relative? Because the variation in the observed values depends on the expected value. Why squared? Well, squaring makes the differences positive, so they don't cancel each other when we add them up. But other than that, there is no special reason to choose the exponent 2. The absolute value would also be a reasonable choice.

For the observed frequencies, the chi-squared statistic is 20.6. By itself, this number doesn't mean anything. We have to compare it to results from the null hypothesis.

Here is the code that generates simulated data:

def SimulateRolls(sides, num_rolls):
"""Generates a Hist of simulated die rolls.

Args:
sides: number of sides on the die
num_rolls: number of times to rolls

Returns:
Hist object
"""
hist = Pmf.Hist()
for i in range(num_rolls):
roll = random.randint(1, sides)
hist.Incr(roll)
return hist

And here is the code that runs the simulations:

count = 0.
num_trials = 1000
num_rolls = 60
threshold = ChiSquared(expected, observed)

for _ in range(num_trials):
simulated = SimulateRolls(sides, num_rolls)
chi2 = ChiSquared(expected, simulated)
if chi2 >= threshold:
count += 1

pvalue = count / num_trials
print 'p-value', pvalue

Out of 1000 simulations, only 34 generated a chi-squared value greater than 20.6, so the estimated p-value is 3.4%. I would characterize that as borderline significant. Based on this result, I think the die is probably crooked, but would I have the guy whacked? No, I would not. Maybe just roughed up a little.

For this problem, I could have computed the sample distribution of the test statistic analytically and computed the p-value quickly and exactly. If I had to run this test many times for large datasets, computational efficiency might be important. But usually it's not.

And accuracy isn't very important either. Remember that the test statistic is arbitrary, and the null hypothesis often involves arbitrary choices, too. There is no point in computing an exact answer to an arbitrary question.

For most problems, we only care about the order of magnitude: if the p-value is smaller that 1/100, the effect is likely to be real; if it is greater than 1/10, probably not. If you think there is a difference between a 4.8% (significant!) and 5.2% (not significant!), you are taking it too seriously.

So the advantages of analysis are mostly irrelevant, but the disadvantages are not:

1) Analysis often dictates the test statistic; simulation lets you choose whatever test statistic is most appropriate. For example, if someone is cheating at craps, they will load the die to increase the frequency of 3 and/or 4, not 1 and/or 6. So in the casino problem the results are suspicious not just because one of the frequencies is high, but also because the frequent value is 3. We could construct a test statistic that captures this domain knowledge (and the resulting p-value would be lower).

2) Analytic methods are inflexible. If you have issues like censored data, non-independence, and long-tailed distributions, you won't find an off-the-shelf test; and unless you are a mathematical statistician, you won't be able to make one. With simulation, these kinds of issues are easy.

3) When people think of analytic methods as black boxes, they often fixate on finding the right test and figuring out how to apply it, instead of thinking carefully about the problem.

In summary, don't worry about finding the "right" test. There's no such thing. Put your effort into choosing a test statistic (one that reflects the effect you are testing) and modeling the null hypothesis. Then simulate it, and count!

Added June 24, 2001: For some additional examples, see this followup post.

Statistics Workshop

2011-05-17T13:03:00.000-07:00

I am leading a statistics workshop based on Think Stats tomorrow. If you want to follow along at home, my notes and the exercises are here.

Think Stats will be published by O'Reilly in June

2011-05-04T17:11:00.000-07:00

No statistics today, just news. I signed a contract with O'Reilly to publish Think Stats. I turn in the manuscript next week, and it should be out in June. So if you find any errors, now is the time to tell me!

Here's my new author page at O'Reilly.

Some people have asked if I get to choose the animal on the cover (which is a signature of O'Reilly cover design). The short answer is no.

First, I'm actually not sure whether Think Stats will be an "animal book" or part of another O'Reilly series. We haven't talked about it.

Also, here's what the O'Reilly author materials have to say on the topic:

Cover Design

The purpose of a book cover is to get a potential reader to pick up the book, and to persuade a bookstore to display it.

We're confident that we have the most striking and effective covers in technical publishing. Despite our relatively small size, it's not unusual to see window displays of our books at technical bookstores.

Our covers are all designed by Edie Freedman. She is open to suggestions, but has the final say on all cover designs. Here's what Edie has to say about how she designs the animal covers for Nutshell Handbooks:

I ask the authors to supply me with a description of the topic of the book. What I am looking for is adjectives that really give me an idea of the "personality" of the topic. Authors are free to make suggestions about animals, but I prefer to deal with adjectives. Once I have the information from the author, I spend some time thinking about it, and then I choose an animal. Sometimes it is based on no more than what the title sounds like. (COFF, for example, sounded like a walrus noise to me.) Sometimes it is very much linked to the name of the book or software. (For example, vi, the "Visual Editor," suggested some beast with huge eyes). And it is always subject to what sort of artwork I can find.
It seems to work best this way. We've tried doing animals that the authors want, and have found that it works better if I do the selection, and then submit it for approval. I don't think anyone has been too disappointed with the final choices we've made.

I'll add that if you do want a particular animal, Edie is more inclined to respond to "right brain" reasons than to obvious mental associations. For example, several people on the net suggested the obvious clam for the Perl handbook, but Edie responded instead to Larry Wall's obscure argument for the camel: ugly but serviceable, able to go long distances without a lot of nourishment.

In any event, you'll have a chance to see and comment on Edie's cover design before we commit it to print. If you really hate it, and can't be persuaded to feel otherwise, she'll probably try again (but no promises!).

So here's what I have for "right brain" ideas: the themes of the books are curiosity and analysis. One of the problems I see with conventional approaches to statistics is that it makes people timid -- afraid to do the wrong thing. One of my goals with this approach is to help students be fearless. So I want something curious and fearless. Ideas?

Here's another article about the animal selection process.

And here's the list of animals they've already used.

And here it is:

The line drawing looks better at higher resolution, but you get the idea. Can anyone name that fish?