## 2.7 Two-Way Analysis of Variance

Let us start by creating a copy of program effort and grouping it into categories 0-4, 5-14 and 15+. This time we will make the copying, recoding and labeling in just one command, using `effort_g` for effort in groups.

```. recode effort (0/4=1 "Weak") (5/14=2 "Moderate") ///
>         (15/max=3 "Strong"),    gen(effort_g) label(effort_g)
(20 differences between effort and effort_g)
```

Here's a table showing steeper declines in countries with strong programs, with a smaller difference between weak and moderate:

```. tabulate effort_g, summarize(change)

rECODE of |
effort | Summary of % Change in CBR between
(Program |            1965 and 1975
Effort) |        Mean   Std. Dev.       Freq.
------------+------------------------------------
Weak |           5           4           7
Moderate |   9.3333333    7.393691           6
Strong |   27.857143   6.3358391           7
------------+------------------------------------
Total |        14.3   11.810343          20
```

Let us create dumy variables for moderate and strong programs.

```. gen effort_mod = effort_g==2 // or effort >=5 & effort < 15

. gen effort_str = effort_g==3 // or effort >=15 & !missing(effort)
```

The comments show how we could create these dummies directly from the original variable, bypassing the creation of the grouped factor.

We are now ready to fit the additive model with weak programs in low setting countries as the reference cell:

```. regress change setting_med setting_high effort_mod effort_str

Source |       SS       df       MS              Number of obs =      20
-------------+------------------------------           F(  4,    15) =   13.55
Model |  2075.80829     4  518.952073           Prob > F      =  0.0001
Residual |   574.39171    15  38.2927807           R-squared     =  0.7833
Total |      2650.2    19  139.484211           Root MSE      =  6.1881

------------------------------------------------------------------------------
change |      Coef.   Std. Err.      T    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
setting_med |  -1.680829   3.854967    -0.44   0.669    -9.897497    6.535839
setting_high |   2.387565    4.45653     0.54   0.600    -7.111304    11.88643
effort_mod |   3.836269   3.574561     1.07   0.300    -3.782727    11.45527
effort_str |    20.6715   4.339232     4.76   0.000     11.42265    29.92036
_cons |   5.379275    3.10526     1.73   0.104     -1.23943    11.99798
------------------------------------------------------------------------------
```

Compare these estimates with the results in Table 2.15 and the anova with Table 2.16.

Countries with strong family planning programs show steeper declines than countries with weak programs at the same level of social setting, on average 21 percentage points more.

This statement is based on the assumption of additivity, namely that the difference in outcomes across categories of program effort is the same at every level of setting. We will test this assumption below.

#### Factor Variables

Before we do that, let us show how we can get exactly the same results, and admitedly with less effort, using factor variables. We simply use the categorical versions of our predictors with the `i.` prefix to instruct Stata to treat them as factors:

```. regress change i.setting_g i.effort_g

Source |       SS       df       MS              Number of obs =      20
-------------+------------------------------           F(  4,    15) =   13.55
Model |  2075.80829     4  518.952073           Prob > F      =  0.0001
Residual |   574.39171    15  38.2927807           R-squared     =  0.7833
Total |      2650.2    19  139.484211           Root MSE      =  6.1881

------------------------------------------------------------------------------
change |      Coef.   Std. Err.      T    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
setting_g |
2  |  -1.680829   3.854967    -0.44   0.669    -9.897497    6.535839
3  |   2.387565    4.45653     0.54   0.600    -7.111304    11.88643
|
effort_g |
2  |   3.836269   3.574561     1.07   0.300    -3.782727    11.45527
3  |    20.6715   4.339232     4.76   0.000     11.42265    29.92036
|
_cons |   5.379275    3.10526     1.73   0.104     -1.23943    11.99798
------------------------------------------------------------------------------
```

You just have to remember that setting 2 and 3 are medium and high, and effort 2 and 3 are moderate and strong programs.

#### Fitted Values

Let us reproduce Table 2.17 in the notes, showing fitted means by setting and effort. I will use `predict` to generate fitted values, and then simply tabulate them by the two relevant factors:

```. predict anovafit
(option xb assumed; fitted values)

. label var anovafit "Two-way anova fit"

. tabulate setting_g effort_g , summarize(anovafit) means

Means of Two-way anova fit

Social |   RECODE of effort (Program
Setting |            Effort)
(Grouped) |      Weak   Moderate     Strong |     Total
-----------+---------------------------------+----------
Low | 5.3792748  9.2155437          . | 7.5714285
Medium | 3.6984456  7.5347152  24.369947 | 8.5999999
High | 7.7668395  11.603108  28.438341 | 23.749999
-----------+---------------------------------+----------
Total | 5.0000001  9.3333331  27.857142 |      14.3
```

Can you get the missing cell in the upper right corner? How about the margins?

#### A Two-Factor Interaction

Let us now consider a model with an interaction between social setting and program effort, so differences in fertility decline by effort will vary by setting. To do this "by hand" we need to create four dummy variables. At this point it becomes hard to create reasonable names in 12 characters or less. I'll use one character for each variable and three for each category. An alternative is to use camelCasing, which saves the spaces used by the underscores.

```. gen se_med_mod = setting_med  * effort_mod

. gen se_med_str = setting_med  * effort_str

. gen se_hi_mod  = setting_high * effort_mod

. gen se_hi_str  = setting_high * effort_str

.
. regress change setting_med setting_high effort_mod effort_str ///
>   se_med_mod se_med_str se_hi_mod se_hi_str
note: se_hi_str omitted because of collinearity

Source |       SS       df       MS              Number of obs =      20
-------------+------------------------------           F(  7,    12) =    8.15
Model |     2189.45     7  312.778571           Prob > F      =  0.0009
Residual |      460.75    12  38.3958333           R-squared     =  0.8261
Total |      2650.2    19  139.484211           Root MSE      =  6.1964

------------------------------------------------------------------------------
change |      Coef.   Std. Err.      T    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
setting_med |   3.333333    5.05937     0.66   0.522    -7.690086    14.35675
setting_high |   6.333333   7.155029     0.89   0.393    -9.256136     21.9228
effort_mod |   8.583333   4.732607     1.81   0.095    -1.728132     18.8948
effort_str |   19.33333   6.692917     2.89   0.014      4.75072    33.91595
se_med_mod |  -14.58333   8.578579    -1.70   0.115    -33.27445    4.107784
se_med_str |  -.3333333   9.797427    -0.03   0.973    -21.68009    21.01343
se_hi_mod |  -6.583333   9.959379    -0.66   0.521    -28.28296    15.11629
se_hi_str |  (omitted)
_cons |   2.666667   3.577515     0.75   0.470    -5.128068     10.4614
------------------------------------------------------------------------------
```

Oops, Stata dropped a term. Why? Because there are no countries with strong programs in low settings, so we have only eight groups, but are trying to represent their means using nine parameters, which is obviously one too many.

Fortunately this doesn't affect testing. We can use the sum of squares to build the hierarchical anova in Table 2.21, and Stata can test the interaction for us, dropping automatically the redundant term:

```. test se_med_mod se_med_str se_hi_mod se_hi_str

( 1)  se_med_mod = 0
( 2)  se_med_str = 0
( 3)  se_hi_mod = 0
( 4)  o.se_hi_str = 0
Constraint 4 dropped

F(  3,    12) =    0.99
Prob > F =    0.4318
```

We have no evidence that the differences by effort vary with social setting.

This makes the issue of interpreting parameters moot, but it may be worth noting briefly the problems caused by the empty cell. As things stand, the coefficient of moderate effort compares moderate with weak at low setting, but the coefficient of strong effort compares strong with weak at high setting. (Table 2.20 in the notes may help see this point. When the term for high and strong is dropped, the only difference between weak and strong programs at high setting is the coefficient of strong.)

The parametrization I like best for this problem combines the main effects of effort with the interactions, so we obtain differences between strong and weak, and between moderate and weak programs, at each level of setting. This allows us to omit the difference between strong and weak programs at low setting, which is the one we can't identify. I will not pursue this tack at this point, but will return to the idea when we have to interpret a significant interaction.

#### Factor Interactions

We can fit the same model using factor variables, we just need to use the `i.` prefix to indicate that grouped setting and effort should be treated as factors, not covariates, combined with a hash `#` to indicate interactions.

Stata understands a single or double hash, typed without spaces after a factor specification:

A single hash, as in `i.setting#i.effort`, takes the first cell in the cross-tabulation of setting and effort groups as the baseline and creates eight dummies, one for each of the other cells. In this parametrization each combination of levels is compared directly to the baseline.

A double hash, as in `i.setting##i.effort`, creaes main effects and interactions, and corresponds to the parametrization we used above. Recall that in this case "main" effects are really differences between categories of one factor when the other factor is at the baseline, and the "interactions" are additional differences when the other factor is not at the baseline.

In our example we have the added complication that one of the cells, (1,3) is empty. Stata detects this and drops the last term:

```. regress change i.setting_g##i.effort_g
note: 1b.setting_g#3.effort_g identifies no observations in the sample
note: 3.setting_g#3.effort_g omitted because of collinearity

Source |       SS       df       MS              Number of obs =      20
-------------+------------------------------           F(  7,    12) =    8.15
Model |     2189.45     7  312.778571           Prob > F      =  0.0009
Residual |      460.75    12  38.3958333           R-squared     =  0.8261
Total |      2650.2    19  139.484211           Root MSE      =  6.1964

------------------------------------------------------------------------------
change |      Coef.   Std. Err.      T    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
setting_g |
2  |   3.333333    5.05937     0.66   0.522    -7.690086    14.35675
3  |   6.333333   7.155029     0.89   0.393    -9.256136     21.9228
|
effort_g |
2  |   8.583333   4.732607     1.81   0.095    -1.728132     18.8948
3  |   19.33333   6.692917     2.89   0.014      4.75072    33.91595
|
setting_g#|
effort_g |
1 3  |  (empty)
2 2  |  -14.58333   8.578579    -1.70   0.115    -33.27445    4.107784
2 3  |  -.3333333   9.797427    -0.03   0.973    -21.68009    21.01343
3 2  |  -6.583333   9.959379    -0.66   0.521    -28.28296    15.11629
3 3  |  (omitted)
|
_cons |   2.666667   3.577515     0.75   0.470    -5.128068     10.4614
------------------------------------------------------------------------------
```

Unfortunately we can't control which of the dummies is dropped. All we can do is change the baseline, or combine categories, but that is not particularly helpful.

Continue with 2.8 Analysis of Covariance Models