Fitting a set of autologistic models Abstract, then setup, then the mathematical details, and full-detail question:
ABSTRACT
Usually, autologistic models are fit to a single graph. I would like to fit the same model to many graphs, simultaneously, because I believe the same "spatial" dependence applies to all of them. What is the best way to do this? (I have some ideas.)
SETUP
I have run surveys on the participants in a large number of discussion groups, resulting in one set of surveys for each of the $K$ distinct groups, with the size $n_k$ of each group ranging from 5-10. In the collected data, I have measures of an outcome $Y_i$ for each participant (their satisfaction with the discussion) and a measure of their perceived similarity to the other participants, $a_{ijk}$, with $A_k$ being the full matrix of perceived participant similarities.
I am assuming that my satisfaction is partially a function of my features $X_i$, and partially the satisfaction of the people I think are like me. Intuitively, I'm not sure how to feel, so I look to my perceived peers for clues. I'd like to fit a model to this process, which will tell me how much my own features matter, and how much I really care about my peers' opinions.
MATH
A traditional autologistic model is a Markov Random Field, meant for a SINGLE network. It requires an adjacency matrix $A$, and an out come vector $Y$, and can handle additional covariates on the nodes $X_i$. It assumes that
$$
\log \frac{P(Y_i = 1|Neighbors(i))}{P(Y_i = 0|Neighbors(i))} = \alpha X_i + \beta \sum_{j \in N(i)} Y_j
$$
Or in other words, that the log odds of your happiness depends on you, and the happiness of those adjacent to you. The parameter $\alpha$ controls the influence of your covariates, and $\beta$ the influence of your neighbors. (Also you should center the models, but I won't muddy the waters here.)
A major problem with these models in most cases is that they're too hard to fit exactly, due to an intractable partition function which requires summing over of $2^N$ states. Unfortunately the approximations that are used instead appear to give good estimates only with very large sample sizes, which I don't have.
What I DO have though, is something in between -- many networks, EACH of small sample size (rendering the partition function tractable), but AGGREGATELy of decent size (possibly making the approximations good). 
My idea was to attempt to fit a model by simply adding the log likelihoods together for each of the groups, and maximizing that directly. This would be done by calculating each of the many, tractable, partition functions within that sum, directly, since there are at most 1024 states to sum over for each one. And at the end, I would just feed it all to an optimizer.
QUESTIONS
Is this reasonable?
If not, what might be a better way to do it?
If it is, are there still things I'm going to need to tell the optimizer to make this work properly, as in, an appropriate choice of algorithm? I have a very feeble grasp of black-box optimizers.
 A: Your suggestion is reasonable.  If your groups can be assumed independent, the joint likelihood of respondents in all groups is just the product of the individual group likelihoods.  
In symbols, if you have response vectors $\mathbf{Y}_1,\mathbf{Y}_2,\ldots \mathbf{Y}_k$ for the $k$ groups, you're assuming that each $\mathbf{Y}_i$ is a vector of correlated binary random variables following an autologistic distribution, call it $g_i(\mathbf{Y};\mathbf{\alpha_i},\beta)$.  Here $g_i(\cdot)$ is the autologistic pmf for the $i$th group (each group could have its own graph), $\mathbf{\alpha}_i$ is the vector of individual-specific parameters for group $i$, and $\beta$ is an association parameter held in common between groups. If the groups are independent the joint log-likelihood is then 
$$l(\mathbf{\alpha}_1,\mathbf{\alpha}_2,\ldots,\mathbf{\alpha}_k,\beta) = \sum_{i=1}^k \log g_i(\mathbf{Y}_i;\mathbf{\alpha_i},\beta).$$
The $g_i(\cdot)$ functions are the actual pmfs, containing the partition functions.  For a single parameter setting, computing all the partition functions means evaluating the un-normalized pmf $\sum_{i=1}^k2^{n_i}$ times, which should be possible (as you mentioned) for your case as long as $k$ isn't very large.
You could feed it to any continuous optimizer (e.g. optim or nlm in R), but it would be nice to double-check whether or not the thing is convex first to know whether you need to worry about getting stuck in local optima.   
Finally, some miscellaneous notes.


*

*You could alternatively think of your data (all groups) as one big graph with $k$ disconnected subgraphs.  The likelihood should work out the same.

*I've written the above as if the $\mathbf{\alpha}_1,\ldots\mathbf{\alpha}_k$ are distinct parameter vectors, but it's possible you meant that the individual-specific part of your model includes a single vector $\mathbf{\alpha}$ and covariate vectors $\mathbf{x}_i$.  This is the autologistic regression model.

*I strongly recommend that you use the $\{-1,1\}$ coding, rather than $\{0,1\}$ values for your responses. It may seem a trivial change but it greatly influences the interpretation of the parameters.  The centered model you mentioned does not do a great job at fixing things.  Better just to use the standard model with plus/minus variable coding.  I have a paper under review explaining this point but as of now I'm not aware of any publicly-available reference about it.

