Previous: transform and derive

With `crunch`

, you can harness the power of R to do
computations with your datasets in Crunch that would be difficult or
impossible to accomplish in a graphical user interface.

While the web application certainly supports crosstabbing, you may
want to do aggregations like this in R. Crosstabbing in R with
`crunch`

may allow you to easily do additional computations
on the result, for example.

`crunch`

contains the `crtabs`

(Crunch-tabs)
function, which largely emulates the design of the `xtabs`

function in base R. In essence, you define a formula and provide data in
which to evaluate it. In this case, we’ll be providing a
`CrunchDataset`

.

Like `xtabs`

, `crtabs`

takes a formula and a
data argument. Dimensions of your crosstab go on the right side of the
`~`

. For a univariate table of frequencies by education, we
can do

```
<- crtabs(~ educ, data=ds)
tab1 tab1
```

```
## educ
## No HS High school graduate Some college 2-year 4-year
## 12 71 61 24 57
## Post-grad
## 25
```

Additional dimensions are added with `+`

. For a two-way
table of education and gender,

```
<- crtabs(~ educ + gender, data=ds)
tab2 tab2
```

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

`crtabs`

takes advantage of several Crunch features that
`xtabs`

does not support. First, it respects weight variables
that have been set on the server. This dataset is not currently
weighted

`weight(ds)`

`## NULL`

but we can very easily change that. Let’s use the “weight” variable that already exists in the dataset:

`weight(ds) <- ds$weight`

Now, if we request the same two-way table as before, we’ll get weighted results:

`crtabs(~ educ + gender, data=ds)`

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

If we want unweighted data, that’s easy enough:

`crtabs(~ educ + gender, data=ds, weight=NULL)`

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

As with any `array`

data type, we can compute margin
tables, and the `prop.table`

function in R provides a
convenient way for sweeping a table by a margin. These work on the
output of `crtabs`

, too:

`prop.table(tab1)`

```
## educ
## No HS High school graduate Some college 2-year 4-year
## 0.048 0.284 0.244 0.096 0.228
## Post-grad
## 0.100
```

For column proportions, specify margin=2 (by rows, margin=1):

`prop.table(tab2, 2)`

```
## gender
## educ Male Female
## No HS 0.05405405 0.04316547
## High school graduate 0.29729730 0.27338129
## Some college 0.23423423 0.25179856
## 2-year 0.05405405 0.12949640
## 4-year 0.22522523 0.23021583
## Post-grad 0.13513514 0.07194245
```

Let’s make that more readable:

`round(100*prop.table(tab2, 2))`

```
## gender
## educ Male Female
## No HS 5 4
## High school graduate 30 27
## Some college 23 25
## 2-year 5 13
## 4-year 23 23
## Post-grad 14 7
```

`crtabs`

also comfortably handles the more complex data
types that Crunch supports, including categorical array and multiple
response variables. In the array
variables vignette, we created a categorical array, “imiss”, for
“Important issues”. We can crosstab with arrays just as we do
non-arrays.

```
<- crtabs(~ imiss + gender, data=ds)
tab3 tab3
```

```
## , , gender = Male
##
## imiss
## imiss Very Important Somewhat Important Not very Important Unimportant
## Abortion 32 31 26 22
## Education 49 44 12 6
## Gay rights 18 25 21 46
## Health care 79 23 5 4
## Immigration 51 38 14 8
## Medicare 58 35 12 6
## Social security 67 33 5 6
## Taxes 75 29 5 1
## Terrorism 44 38 17 11
## The budget deficit 60 26 16 9
## The economy 93 14 4 0
## The environment 35 42 20 14
## The war in Afghanistan 24 51 25 10
##
## , , gender = Female
##
## imiss
## imiss Very Important Somewhat Important Not very Important Unimportant
## Abortion 63 42 24 9
## Education 87 41 7 4
## Gay rights 41 36 25 37
## Health care 102 30 4 2
## Immigration 56 53 20 9
## Medicare 80 41 15 2
## Social security 82 41 14 2
## Taxes 85 35 15 4
## Terrorism 68 44 18 9
## The budget deficit 65 50 17 7
## The economy 103 30 4 1
## The environment 64 45 16 11
## The war in Afghanistan 49 57 20 13
```

Note that even though we specified two variables in our formula, because “imiss” itself is two dimensional, our result is a three-dimensional array.

To illustrate working with multiple response variables, let’s convert “imiss” to multiple response, selecting its positive categories as indicating selection:

`$imiss <- dichotomize(ds$imiss, c("Very Important", "Somewhat Important")) ds`

Now, when we crosstab it, we’ll get a two-dimensional table because multiple response variables present a one-dimensional interface:

```
<- crtabs(~ imiss + gender, data=ds)
tab3mr tab3mr
```

```
## gender
## imiss Male Female
## Abortion 58.64993 97.61001
## Education 93.98404 115.08118
## Gay rights 39.50817 68.21617
## Health care 115.13064 120.27883
## Immigration 86.65671 96.90796
## Medicare 107.80704 110.61446
## Social security 109.95809 109.20223
## Taxes 106.60644 108.92159
## Terrorism 88.26523 101.57750
## The budget deficit 89.38060 104.02191
## The economy 108.38178 119.16976
## The environment 90.32808 97.05483
## The war in Afghanistan 73.62479 96.68801
```

It’s worth noting here that the result of `crtabs`

isn’t
an `array`

object but a `CrunchCube`

object.

`class(tab3mr)`

```
## [1] "CrunchCube"
## attr(,"package")
## [1] "crunch"
```

This allows us to do the appropriate calculations on arrays and
multiple response variables when `prop.table`

is called. To
compute a margin table over a multiple response variable, summing along
the dimension would give an incorrect value because the responses in a
multiple response are not mutually exclusive–they can’t be assumed to
sum to 100 percent. However, the `margin.table`

method on
`CrunchCubes`

can compute the correct margin, so
`prop.table`

gives correct proportions:

`round(100*prop.table(tab3mr, 2))`

```
## gender
## imiss Male Female
## Abortion 48 77
## Education 77 90
## Gay rights 33 53
## Health care 95 94
## Immigration 71 77
## Medicare 89 87
## Social security 90 85
## Taxes 88 85
## Terrorism 74 79
## The budget deficit 74 81
## The economy 89 93
## The environment 74 78
## The war in Afghanistan 61 75
```

Finally, just as we saw in the array variables vignette, we can grab individual subvariables and crosstab with them:

`crtabs(~ imiss$imiss_f + gender, data=ds)`

```
## gender
## imiss_f Male Female
## Very Important 44 68
## Somewhat Important 38 44
## Not very Important 17 18
## Unimportant 11 9
```

It’s worth noting that we can extend the crosstabbing to higher dimensions, just by adding more terms on the right-hand side of the formula:

`round(crtabs(~ imiss + educ + gender, data=ds))`

```
## , , gender = Male
##
## educ
## imiss No HS High school graduate Some college 2-year 4-year Post-grad
## Abortion 3 19 18 1 14 8
## Education 6 23 22 6 23 13
## Gay rights 1 11 12 1 8 10
## Health care 5 33 22 6 22 14
## Immigration 4 26 20 6 20 13
## Medicare 5 29 21 5 21 12
## Social security 5 31 22 6 21 15
## Taxes 6 30 25 6 23 14
## Terrorism 5 24 19 6 16 12
## The budget deficit 6 25 20 5 18 12
## The economy 6 31 25 6 24 15
## The environment 5 22 16 4 17 13
## The war in Afghanistan 5 17 20 3 20 10
##
## , , gender = Female
##
## educ
## imiss No HS High school graduate Some college 2-year 4-year Post-grad
## Abortion 4 26 27 12 27 9
## Education 6 31 32 17 32 10
## Gay rights 3 15 23 9 21 6
## Health care 6 35 33 18 30 10
## Immigration 5 27 30 13 24 10
## Medicare 6 32 30 18 27 8
## Social security 6 33 33 16 27 8
## Taxes 6 30 29 15 31 9
## Terrorism 4 32 27 16 25 8
## The budget deficit 5 32 28 16 26 8
## The economy 6 34 33 18 32 10
## The environment 5 25 24 16 30 9
## The war in Afghanistan 5 26 31 15 22 7
```

`crtabs`

can also compute quantities other than counts.
Using the left-hand side of the formula, we can specify other
aggregations to put in the cells of the table. For example, in the deriving variables vignette, we created an “age”
variable. We can easily compute the average age by gender and
education:

`crtabs(mean(age) ~ educ + gender, data=ds)`

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

Other supported aggregations include `min`

,
`max`

, `sd`

, and `sum`

. For the minimum
age by gender and education,

`crtabs(min(age) ~ educ + gender, data=ds)`

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

We can get unconditional (univariate) statistics by making the
right-hand side of your formula be just the number `1`

:

`crtabs(min(age) ~ 1, data=ds)`

`## [1] 250`

Numeric aggregation functions also work with categorical variables that have numeric values defined for their categories; this is the reason why numeric values for categories are defined, in fact. In the variables vignette, we worked with the “On the right track” question and set some numeric values:

`categories(ds$track)`

```
## id name value missing
## 1 1 Generally headed in the right direction 1 FALSE
## 2 3 Not sure 0 TRUE
## 3 2 Wrong track -1 FALSE
## 4 -1 No Data NA TRUE
```

We can use these numeric values to compute an “on the right track index” by averaging them. If the index is greater than zero, more people thing things are going well, and if it is negative, more respondents are pessimistic.

`round(crtabs(mean(track) ~ educ + gender, data=ds), 2)`

```
## gender
## educ Male Female
## No HS 6 6
## High school graduate 33 38
## Some college 26 35
## 2-year 6 18
## 4-year 25 32
## Post-grad 15 10
```

Looks like most people surveyed thought that the U.S. is on the wrong track, but that pessimism is less pronounced for women with higher levels of education.

We can also specify a subset of `ds`

to analyze, just as
if it were a data.frame. Let’s do the same calculation for Democrats
only:

`round(crtabs(mean(track) ~ educ + gender, data=ds[ds$pid3 == "Democrat",]), 2)`

```
## gender
## educ Male Female
## No HS 0 3
## High school graduate 7 17
## Some college 8 16
## 2-year 1 6
## 4-year 9 19
## Post-grad 5 5
```

Not surprisingly, Democrats were less pessimistic about the direction of the country than the general population.

A few final observations about `crtabs`

. First, all of
these calculations have been weighted by the weight variable we set
above. We set it and could then forget about it–we didn’t have to litter
all of our expressions with `ds$weight`

and extra arithmetic
to do the weighting. Crunch handles this for us.

Second, none of these aggregations required pulling case-level data
to your computer. `crtabs`

sends Crunch expressions to the
server and receives in return an `n`

-D array of results. The
only computations happening locally are the margin tables and sweeping
in `prop.table`

, computing on the aggregate results. Your
computer would work exactly as hard with this example dataset of 1000
rows as it would with a dataset of 100 million rows.

Any statistical modeling function that takes a `data`

argument should happily accept a `CrunchDataset`

and just do
the right thing–no extra effort or thought required.

Let’s fit a basic Ordinary Least Squares (OLS) model. In our dataset, we have a few questions about Edward Snowden, such as:

`$snowdenleakapp ds`

```
## Approval of Snowden's Leak (categorical)
##
## Count
## Strongly disapprove 86.00644
## Not sure 55.82542
## Somewhat approve 43.19385
## Somewhat disapprove 40.26587
## Strongly approve 24.70842
```

We can use `lm`

to fit our model. Let’s explore the
relationship between approval of Snowden’s leak and respondents’
interest in current events, party identification, gender, and age.

```
<- lm(I(snowdenleakapp == "Strongly approve") ~ newsint2 + pid3 + gender + age,
ols1 data=ds)
summary(ols1)
```

```
##
## Call:
## lm(formula = I(snowdenleakapp == "Strongly approve") ~ newsint2 +
## pid3 + gender + age, data = ds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33261 -0.15753 -0.12930 -0.06006 0.94924
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.133723 0.093700 1.427 0.155
## newsint2Some of the time 0.016157 0.050750 0.318 0.750
## newsint2Only now and then -0.007068 0.074048 -0.095 0.924
## newsint2Hardly at all -0.058744 0.091630 -0.641 0.522
## pid3Republican -0.088822 0.058002 -1.531 0.127
## pid3Independent 0.022139 0.052650 0.420 0.675
## pid3Other 0.185892 0.198541 0.936 0.350
## pid3Not sure -0.049288 0.103442 -0.476 0.634
## genderFemale -0.019411 0.045856 -0.423 0.672
## age 0.000361 0.001651 0.219 0.827
##
## Residual standard error: 0.3366 on 240 degrees of freedom
## Multiple R-squared: 0.02533, Adjusted R-squared: -0.01122
## F-statistic: 0.6931 on 9 and 240 DF, p-value: 0.7149
```

Looks like partisanship is weakly associated with approval of the NSA
leak, but overall the model isn’t a great fit, given our data. (For what
it’s worth, we’re working on a randomly drawn subset of the survey so
that the size of data included with package is small. Results are more
meaningful with the full dataset.) Nevertheless, this example
illustrates how straightforward it is to do statistical analysis with
data in Crunch. Even though your dataset lives on the server, you can
think of it like a local `data.frame`

. Note, for example,
that our categorical variables (News Interest, Party ID, and Gender)
expand their categories out as dichotomous indicators, just as if they
were `factor`

variables in a `data.frame`

.

Given that we’re estimating a model with a dichotomous dependent
variable, perhaps a logit would be more appropriate than a strict linear
predictor. We can use `glm`

instead:

```
<- glm(I(snowdenleakapp == "Strongly approve") ~ newsint2 + pid3 + gender + age,
logit1 family=binomial(link="logit"), data=ds)
summary(logit1)
```

```
##
## Call:
## glm(formula = I(snowdenleakapp == "Strongly approve") ~ newsint2 +
## pid3 + gender + age, family = binomial(link = "logit"), data = ds)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8982 -0.5873 -0.5249 -0.3413 2.4261
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.883494 0.847876 -2.221 0.0263 *
## newsint2Some of the time 0.141215 0.439458 0.321 0.7480
## newsint2Only now and then -0.075103 0.700499 -0.107 0.9146
## newsint2Hardly at all -0.734607 1.084190 -0.678 0.4980
## pid3Republican -1.077062 0.673224 -1.600 0.1096
## pid3Independent 0.156299 0.439101 0.356 0.7219
## pid3Other 1.056509 1.273998 0.829 0.4069
## pid3Not sure -0.535621 1.119244 -0.479 0.6323
## genderFemale -0.176401 0.410364 -0.430 0.6673
## age 0.003547 0.015024 0.236 0.8134
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 191.28 on 249 degrees of freedom
## Residual deviance: 184.52 on 240 degrees of freedom
## AIC: 204.52
##
## Number of Fisher Scoring iterations: 5
```

As before, not a particularly interesting result, but this is just
the beginning of the analysis process. Using `crunch`

, you
can keep exploring the data and perhaps find a better fit.

Unlike the previous examples, these modeling functions do have to
pull columns of data from the server to your local machine. However,
only the columns of data you reference in your formula are copied, and
if you specify a subset of the dataset to regress on (as we did above
with `crtabs`

when we looked at just Democrats), only those
rows are retrieved. This helps minimize the time spent shipping data
across the network. Moreover, because of the `crunch`

package’s query cache, subsequent models that incorporate any of those
variables will not have to go to the server to get them.