# Interpreting GLS in the context of modeling a multivariate distribution

I have some paired observations $$(x,y)$$ which form some type of distribution. In an attempt to simplify things, I'm trying to fit a distribution of $$Y=\beta_0 + \beta_1*X + Error$$ (linear regression) and fit the parameters by maximum likelihood. My idea is that $$P(X)*P(Y|X)=P(Y)*P(X|Y)=P(X,Y)$$ where I use generalized least squares to calculate $$P(Y|X)$$. I used generalized least squares because the error has a special correlation structure. In my model, the error structure can be one of two different models, the distribution of $$X$$ can be one of two distributions, and the distribution of $$Y$$ can be one of two distributions. So I know one of the 8 combinations of distributions is correct for the data. However, whatever distribution of error structure I use in the GLS error matrix seems to force me to also model $$Y$$ and $$X$$ as having the same distribution also. What am I doing wrong?

If it helps, the exact details are that I'm modeling evolution. Species have trait values based on their length of shared evolutionary history. Because woodpeckers only evolved 20,000,000 years ago they have pretty similar wing structure, compared to a woodpecker versus a hummingbird. So you can see that more closely related species have more closely related traits. However, some traits have opposite patterns. For example, woodpeckers which live together actually have different diets because if they had the same diet whichever was better at exploiting it would drive the other extinct. So for that trait species which are more closely related have less similarity. In theory, depending on our hypotheses, traits can have any level of similarity. I want to run a regression where mouth size = diet + body size, but I only have data for body size and mouth size, so diet is modeled as error. I want to test different models where the three traits have different levels of similarity. But every time I specify a covariance structure for the error it seems to force the dependent trait to also be modeled with that same covariance structure which means I can't test hypotheses where it has a different level of similarity. To drive the point home, notice that for a GLS regression of Y and X

Likelihood(regression of Y on X with GLS error structure)+Likelihood(X modeled as normal distribution) does NOT equal Likelihood(regression of X on Y with same GLS error structure)+Likelihood(Y modeled as normal distribution).

To me this proves the dependent variable cannot be modeled as coming from a normal distribution, or else the likelihood wouldn't change when we switch the variable places. What is going on?!