Specifying a covariance structure: pros and cons What are the benefits of specifying a covariance structure in a GLM (rather than treating  all off-diagonal entries in the covariance matrix as zero)? Aside from reflecting what one knows of the data, does it


*

*improve goodness of fit?

*improve predictive accuracy on held-out data?

*allow us to estimate the extent of covariance?


What are the costs of imposing a covariance structure? Does it


*

*add computational complications for estimation algorithms?

*increase the number of estimated parameters, also increasing AIC, BIC, DIC?


Is it possible to determine the right covariance structure empirically, or is this something that depends on your knowledge of the data-generative process?
Any costs / benefits I didn't mention?
 A: Here's another incomplete answer that isn't even directly about GLM...In my very limited experience with structural equation modeling (SEM), I've picked up a couple ideas that I hope might add something to the discussion. Please bear in mind throughout that I'm speaking from (limited) experience with SEM, not GLM per se, and I'm fairly ignorant of whether and where this distinction might become important. I'm more of a stats user than a statistician, so I'm also not sure that these ideas will apply to all or even most data; I've only found that they have applied to most of my own. 
First, I'd echo @StephanKolassa's emphasis on the importance of modeling what you already know. You acknowledge this as an aside, but I think the benefits that you're asking about are benefits of modeling what you know. As such, they meaningfully reflect that your resultant model possesses the information about the covariance structure that you've added.  
In SEM, I have found (through limited experience, not through theoretical study):
Benefits


*

*Modeling the covariance structure improves goodness of fit (GoF) if the covariance is much stronger than its standard error (i.e., if the symmetric pathway is significant). This means you usually won't improve GoF by modeling near-zero correlations, and multicollinearity can cause problems for GoF because it inflates standard errors.

*Haven't tried holding out data to predict yet, but my intuition is that fixing the covariances to zero in your model is analogous to predicting a DV by combining a set of separate, single-IV, linear regression equations. Unlike this approach, multiple regression accounts for covariance in the IVs when producing a model of equations to predict the DV. This certainly improves interpretability by separating direct effects from indirect effects that occur entirely within the included set of IVs. Honestly, I'm not sure whether this necessarily improves prediction of the DV though. Being a stats user and not a statistician, I threw together the following simulation testing function to give an incomplete answer (apparently, "Yes, predictive accuracy improves when the model incorporates IV covariance") in this hopefully analogous case...
simtestit=function(Sample.Size=100,Iterations=1000,IV.r=.3,DV.x.r=.4,DV.z.r=.4) {
require(psych); output=matrix(NA,nrow=Iterations,ncol=6); for(i in 1:Iterations) {
x=rnorm(Sample.Size); z=rnorm(Sample.Size)+x*IV.r
y=rnorm(Sample.Size)+x*DV.x.r+z*DV.z.r
y.predicted=x*lm(y~x+z)$coefficients[2]+z*lm(y~x+z)$coefficients[3]
bizarro.y.predicted=x*lm(y~x)$coefficients[2]+z*lm(y~z)$coefficients[2]
output[i,]=c(cor(y.predicted,y)^2,cor(bizarro.y.predicted,y)^2,
cor(y.predicted,y)^2>cor(bizarro.y.predicted,y)^2,cor(x,z),cor(x,y),cor(y,z))}
list(output=output,percent.of.predictions.improved=100*sum(output[,3])/Iterations,
mean.improvement=fisherz2r(mean(fisherz(output[,1])-fisherz(output[,2]))))}

# Wrapping the function in str( ) gives you the gist without filling your whole screen
str(simtestit())

This function generates random samples ($N =$ Iterations, $n$ = Sample.Size) from three normally distributed variables: z $=$ x $+$ random noise, and y $=$ x $+$ z $+$ random noise. The user can influence their correlations somewhat by overriding the defaults for the last three arguments, but the random noise affects the sample correlations too, so this simulates the way sampling error affects estimates of true correlation parameters. The function computes predictions of y based on regression coefficients for x and z derived from:
($1$) multiple regression (y.predicted), and...
($2$) two separate, bivariate linear regressions (bizarro.y.predicted). 
The output matrix contains Iterations rows and six columns: the $R^2$s of $1$ and $2$, a true-false test of whether $1 > 2$, and the bivariate $r$s for the three unique combinations of x, y, & z. This function produces a three-element list, the first of which is the output matrix. By default, this is 1,000 rows long, so I recommend wrapping simtestit() in the str( ) function or removing this element from the list in the function itself unless you're interested in the individual sample stats for some reason. The the percentage of iterations in which $R^2$ was improved by using ($1$) multiple regression to account for the covariance of the IVs, and the mean of these improvements across the iterations (in the scale of $r$, using a Fisher transformation via the psych package).
The function defaults to a short sim test of fairly typical circumstances for a maximally basic multiple regression. It permits the user to change individual sample sizes and variable correlations to suit the study and prior theories of relationship strength. I haven't tested all possible settings, but every time I've run the function, 100% of the iterations have produced higher $R^2$ with multiple regression. The mean improvement in $R^2$ seems to be greater when the covariance of the IVs (which can be manipulated incompletely by entering an argument for IV.r) is larger. Since you're probably more familiar with your GLM function than I am (which is not at all), you could probably change this function or use the basic idea to compare GLM predictions across however many IVs you want without too much trouble. Assuming that would (or does) turn out the same way, it would seem that the basic answer to your second question is probably yes, but how much depends on how strongly the IVs covary. Differences in sampling error between the held-out data and the data used to fit the model could overwhelm the improvement in its predictive accuracy within the latter dataset, because again, the improvement seems to be small unless IV correlations are strong (at least, in the maximally basic case with only two IVs).

*Specifying a free path for covariance between IVs in the model means asking the model fitting function to estimate this pathway's coefficient, which represents the extent of covariance between IVs. If your GLM function allows you to specify a model in which the the covariance between the IVs is freely estimated rather than fixed to zero, then your problem is a hopefully simple matter of figuring out how to do this and how to get your function to output that estimate.  If your function estimates IV covariances by default, your problem simplifies further to just the latter matter (as is the case with lm( )).
Costs


*

*Yes, freely estimating covariance between IVs means the model fitting algorithm has to do some work to estimate that pathway's coefficient.  Not specifying that pathway in the model usually means fixing the coefficient to zero, which means the model fitting algorithm doesn't need to estimate the coefficient.  Estimating additional covariance parameters means the overall model will require more time to fit.  In models that already take a long time to estimate, the extra time can be substantial, especially if you have a lot of IVs.

*Yes, a freely-estimated covariance structure implies parameter estimates. Populations have covariance parameters, so if you're estimating population covariances, you're estimating parameters. However, if your model fits much better because you're choosing to estimate a non-trivial correlation rather than fixing it to zero, you can probably expect the Akaike and Bayesian information criteria to improve, just like other criteria that incorporate GoF. I'm not familiar with the deviance information criterion (the DIC to which you're referring, right?), but judging from its Wikipedia page, it also seems to incorporate GoF and a penalty for model complexity. 
Therefore the GoF should just need to improve proportionally more than the model's complexity increases to improve the DIC. If this doesn't happen overall, criteria like these that penalize for model complexity will worsen as you estimate more IV covariances. This could be a problem if, for instance, your IVs don't correlate, but the covariance structure is freely estimated anyway because you think the IVs might correlate, or because that's your function's default setting. If you have prior theoretical reasons to assume a correlation is zero and you don't want your model to test this assumption, this is one case where you might be justified in fixing the path to zero. If your prior theory is approximately right, indices that penalize for model complexity will improve if you fix pathways to your prior theory instead of having the model fitting algorithm estimate them freely.
Dunno which function you're working with, but once again, I'm sure I'm unfamiliar with it, so I'm sure this answer could be improved, especially my answer to the second benefit question (for one thing, a mathematical proof of what I'm answering by simulation about multiple regression is probably available somewhere out there). I'm not even familiar with GLM in general (assuming you do mean generalized, not general linear modeling, as the tag suggests), so I hope someone will comment on or edit this answer if the distinctions from SEM invalidate my answers to your questions at all. 
Nonetheless, it seems we've been waiting ten months for the gurus to speak up, so if this doesn't get them to do it, it'll just have to do by itself, I suppose. Let me know if you have a particular GLM function in mind that you'd like me to mess with in R though. I may be able to figure out how to answer #3 more directly for your application if you can specify a GLM function of interest in R. I'm no expert with simulation testing either, but I think your other four questions could be sim tested (more directly) too.
A: Basically, you must specify a covariance structure in GLM. If by "assuming no covariance", you mean "all off-diagonal entries in the covariance matrix are zero", then all you did was assume one very specific covariance structure. (You could be even more specific, e.g., by assuming that all variances are equal.)
This is really a variation on "I don't subscribe to any philosophy; I'm a pragmatist." - "You just described the philosophy you subscribe to."
As such, I would say that the advantage of thinking about the covariance structure is the chance of using a model that is more appropriate to your data. Just as you should include known functional relationships for the expected value (or the mean) of your observations, you should account for any structure you know in the covariance.
And of course, the "disadvantage" is that you need to actually think about all this. Much easier to just use your software's default setting. But this is kind of like always driving in the first gear because your car was in first gear when you bought it and understanding the gear shift takes effort. Not recommended.
