Grouping probabilities of full model vs fitting model on subset of variables
A very thought provoking question. 

Surprisingly, second approach (subset) is better in theory. The first one (full) is an unbiased estimator of second one (subset). That is, expectation of average of probabilities in full model is equal to the probability in subset model.

I said "in theory", because we assume that full and subset models estimate the true P(event | full) and P(event | subset) exactly respectively (zero error, perfect generalization, etc.). In practice, choice of learning algorithm and training set affects how accurate they can fulfill this assumption. A model could be Naive Bayes, Logistic Regression, Neural Network with a SoftMax layer, etc.

Proof

Let 


$F_f$ denote the full features, 
$F_s$ denote the subset features (that we group on), 
$F_d=F_f-F_s$ denote the features that we average out, 
$\boldsymbol{A}=\{x|F_s(x)=a\}$ denote the set of all possible instances with subset features equal to $a$, and 
$A \subset \boldsymbol{A}$ denote the observed subset of $\boldsymbol{A}$ in training set. 


Here is the average of probabilities in full model (first approach):
$$\begin{align*}
P_f(e|F_s=a)&amp;=\frac{1}{N_A}\sum_{x\in A}P(e|F_d(x),a)
\end{align*}$$

Expectation of above average is:
$$\begin{align*}
E[P_f(e|F_s=a)]&amp;=\frac{1}{N_A}\sum_{x\in A}E[P(e|F_d(x),a)]\\
&amp;=\frac{1}{N_A} N_A E[P(e|F_f(X))|X \in \boldsymbol{A}]\\
&amp;=\sum_{b}P(e|F_d=b,F_s=a)P(F_d=b|F_s=a)\\
&amp;=\sum_{b}P(e,F_d=b|F_s=a)\\
&amp;=P(e|F_s=a)
\end{align*}$$

where $P(e|F_s=a)$ is the output of second approach.

Example 

To verify, here is an example:

$e=\text{has_job}$

$F_f=\{\text{sex}, \text{has_degree}\}$

$F_s=\{\text{sex}\}$

$F_d=\{\text{has_degree}\}$

Full model sees:

x     sex  has_degree  has_job
1     f    1           0
2     f    1           1
3     f    0           0

4     m    0           0
5     m    0           0
6     m    1           1


And builds

sex  has_degree  P(has_job=1|sex, has_degree)
f    0           0
f    1           0.5
m    0           0
m    1           1.0


And answers $P(\text{has_job}=1|\text{sex}=f)$ with 

$$\frac{1}{3}\sum_{x \in \{1, 2, 3\}}P(\text{has_job}=1|x)=(0.5+0.5+0)/3=0.33,$$

and answers $P(\text{has_job}=1|\text{sex}=m)$ with 

$$\frac{1}{3}\sum_{x \in \{4, 5, 6\}}P(\text{has_job}=1|x)=(0+0+1.0)/3=0.33$$

On the other hand, subset model sees:

x     sex  has_job
1     f    0
2     f    1
3     f    0

4     m    0
5     m    0
6     m    1


And builds

sex  P(has_job=1|sex)
f    0.33
m    0.33


And answers $P(\text{has_job}=1|\text{sex}=f)$ with 0.33, and $P(\text{has_job}=1|\text{sex}=m)$ with 0.33.