# Efficiently processing a large MxNx2 logistic regression, only interactions matter

I'm working with a large 3-way contingency table (roughly $180 \times 40 \times 2$) — both independent variables are categorical and the response is binary. One independent variable (X) represents a classification of experimental test objects into types, and the other (Y) a set of test conditions. The experiment is only balanced for Y - that is, each test object was tested in all Y conditions, but each test object is only one X type, and the number of objects assigned to each type varies by three orders of magnitude (40­–69600). Response to the test is either "normal" or "anomalous", and, critically, I know a priori that only interactions between X and Y are interesting. I know it's necessary to model the main effects of X and Y as well, but (in the context of the larger experiment) they are noise factors. So the thing that I care about is identifying all pairs of levels of X and Y such that the interaction is significant (given some threshold).

In R-ese, the natural model is

glm(cbind(anom, total-anom) ~ X * Y, family="binomial", data=apt.d)


where apt.d (a $10\times10\times2$ downsample) is at the end of this question. Now, I have two problems:

1. summary(fit object) tells me that this model leaves me with zero residual degrees of freedom and a residual deviance of $4 \times 10^{-10}$, which makes me worry about overfitting, but I don't see how I could structure the model any other way.

2. Actually running the fit on the full data frame is painfully slow -- nearly an hour just to generate the fit object. And then it says that the fit didn't converge and "fitted probabilities numerically 0 or 1 occurred". Is there an alternative implementation of logistic regression in R that would be more appropriate for this size of data set and/or this sort of unbalanced data with widely varying levels of anomalous response?

Downsampled data:

apt.d <- structure(list(X = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L), .Label = c("A",
"B", "C", "D", "E", "F", "G", "H", "I", "J"), class = "factor"),
Y = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L,
5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,
9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L,
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 1L, 2L, 3L, 4L,
5L, 6L, 7L, 8L, 9L, 10L), .Label = c("A", "B", "C", "D",
"E", "F", "G", "H", "I", "J"), class = "factor"), total = c(765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L, 765L,
1530L, 894L, 49L, 77L, 97L, 779L, 1111L, 445L, 45L), anom = c(148L,
235L, 150L, 3L, 0L, 12L, 113L, 148L, 74L, 1L, 121L, 290L,
116L, 1L, 3L, 15L, 276L, 220L, 44L, 2L, 47L, 66L, 23L, 1L,
0L, 3L, 50L, 28L, 20L, 4L, 51L, 42L, 64L, 13L, 0L, 3L, 33L,
67L, 42L, 0L, 91L, 148L, 128L, 0L, 0L, 10L, 69L, 106L, 53L,
4L, 242L, 546L, 307L, 13L, 14L, 25L, 320L, 431L, 106L, 22L,
45L, 37L, 61L, 11L, 0L, 3L, 20L, 28L, 28L, 0L, 143L, 251L,
150L, 16L, 12L, 26L, 118L, 193L, 61L, 3L, 34L, 58L, 29L,
7L, 0L, 3L, 40L, 36L, 33L, 0L, 43L, 32L, 73L, 1L, 1L, 3L,
24L, 30L, 24L, 0L), p.anom = c(0.193464052287582, 0.15359477124183,
0.167785234899329, 0.0612244897959184, 0, 0.123711340206186,
0.145057766367137, 0.133213321332133, 0.166292134831461,
0.0222222222222222, 0.158169934640523, 0.189542483660131,
0.129753914988814, 0.0204081632653061, 0.038961038961039,
0.154639175257732, 0.354300385109114, 0.198019801980198,
0.098876404494382, 0.0444444444444444, 0.061437908496732,
0.0431372549019608, 0.0257270693512304, 0.0204081632653061,
0, 0.0309278350515464, 0.0641848523748395, 0.0252025202520252,
0.0449438202247191, 0.0888888888888889, 0.0666666666666667,
0.0274509803921569, 0.0715883668903803, 0.26530612244898,
0, 0.0309278350515464, 0.0423620025673941, 0.0603060306030603,
0.0943820224719101, 0, 0.118954248366013, 0.0967320261437909,
0.143176733780761, 0, 0, 0.103092783505155, 0.0885750962772786,
0.0954095409540954, 0.119101123595506, 0.0888888888888889,
0.316339869281046, 0.356862745098039, 0.343400447427293,
0.26530612244898, 0.181818181818182, 0.257731958762887, 0.410783055198973,
0.387938793879388, 0.238202247191011, 0.488888888888889,
0.0588235294117647, 0.0241830065359477, 0.0682326621923937,
0.224489795918367, 0, 0.0309278350515464, 0.0256739409499358,
0.0252025202520252, 0.0629213483146067, 0, 0.186928104575163,
0.164052287581699, 0.167785234899329, 0.326530612244898,
0.155844155844156, 0.268041237113402, 0.151476251604621,
0.173717371737174, 0.137078651685393, 0.0666666666666667,
0.0444444444444444, 0.0379084967320261, 0.0324384787472036,
0.142857142857143, 0, 0.0309278350515464, 0.0513478818998716,
0.0324032403240324, 0.0741573033707865, 0, 0.0562091503267974,
0.0209150326797386, 0.08165548098434, 0.0204081632653061,
0.012987012987013, 0.0309278350515464, 0.030808729139923,
0.027002700270027, 0.0539325842696629, 0)), .Names = c("X",
"Y", "total", "anom", "p.anom"), row.names = c(NA, -100L), class = "data.frame")