Logistic regression will, up to numerical imprecision, give exactly the same fits as the tabulated percentages. Therefore, if your independent variables are factor objects factor1
, etc., and the dependent results (0 and 1) are x
, then you can obtain the effects with an expression like
aggregate(x, list(factor1, <etc>), FUN=mean)
Compare this to
glm(x ~ factor1 * <etc>, family=binomial(link="logit"))
As an example, let's generate some random data:
set.seed(17)
n <- 1000
x <- sample(c(0,1), n, replace=TRUE)
factor1 <- as.factor(floor(2*runif(n)))
factor2 <- as.factor(floor(3*runif(n)))
factor3 <- as.factor(floor(4*runif(n)))
The summary is obtained with
aggregate.results <- aggregate(x, list(factor1, factor2, factor3), FUN=mean)
aggregate.results
Its output includes
Group.1 Group.2 Group.3 x
1 0 0 0 0.5128205
2 1 0 0 0.4210526
3 0 1 0 0.5454545
4 1 1 0 0.6071429
5 0 2 0 0.4736842
6 1 2 0 0.5000000
...
24 1 2 3 0.5227273
For future reference, the estimate for factors at levels (1,2,0) in row 6 of the output is 0.5.
The logistic regression gives up its coefficients this way:
model <- glm(x ~ factor1 * factor2 * factor3, family=binomial(link="logit"))
b <- model$coefficients
To use them, we need the logistic function:
logistic <- function(x) 1 / (1 + exp(-x))
To obtain, e.g., the estimate for factors at levels (1,2,0), compute
logistic (b["(Intercept)"] + b["factor11"] + b["factor22"] + b["factor11:factor22"])
(Notice how all interactions must be included in the model and all associated coefficients have to be applied to obtain a correct estimate.)
The output is
(Intercept)
0.5
agreeing with the results of aggregate
. (The "(Intercept)" heading in the output is a vestige of the input and effectively meaningless for this calculation.)
The same information in yet another form appears in the output of table
. E.g., the (lengthy) output of
table(x, factor1, factor2, factor3)
includes this panel:
, , factor2 = 2, factor3 = 0
factor1
x 0 1
0 20 21
1 18 21
The column for factor1
= 1 corresponds to the three factors at levels (1,2,0) and shows that $21/(21+21) = 0.5$ of the values of x
equal $1$, agreeing with what we read out of aggregate
and glm
.
Finally, a combination of factors yielding the highest proportion in the dataset is conveniently obtained from the output of aggregate
:
> aggregate.results[which.max(aggregate.results$x),]
Group.1 Group.2 Group.3 x
4 1 1 0 0.6071429