GLS (nlme) model with interaction and unbalanced data

Question

I have a dataset of soil pH data and two explanatory factors, one with 5 levels (land use)and another one with >100 levels (soil type). Data exploration and comparison of several GLS models fit with nlme suggest that there is an interaction between the two factors. I use GLS instead of lm because there seems to be some heterogeneity in the residuals depending on Soil Type, that I would like to fix with GLS (varIdent). I can test the significance of the interaction with the likelihood-ratio test (ML), but I don't know how to include the interaction term and the vaIdent (by Soil type) with the unbalanced data. Ideally, I would write the model like this:

gls0 <- gls(pH ~ SoilType * LandUse, weights = varIdent(form=~1|SoilType), data=mydata)

But it gives the following error (because the data is unbalanced):

  computed "gls" fit is singular, rank 312

So I create a new factor with the interaction, fit the model, and check that the interaction term is significant:

mydata$ii <- droplevels(interaction(mydata$SoilType, mydata$LandUse))

gls.ii <- gls(pH ~ ii, data = mydata, na.action = na.omit)
gls.m1 <- gls(pH ~ SoilType + LandUse, data=mydata)
anova(update(gls.ii, method = "ML"),
      update(glsm1, method = "ML"))
# or
gls.m1 <- gls(pH ~ SoilType, data=mydata)
anova(update(gls.ii, method = "ML"),
      update(glsm1, method = "ML"))

My first isssue is that I cannot fit the model gls.ii <- gls(pH ~ ii, weights = varIdent(form=~1|ii), data = mydata, na.action = na.omit) for memory issues (Error: cannot allocate vector of size 17.8 Gb). This is why I tested the interaction term without including the error variance term in the models.

But my second question is, is there another way I can fit the model with unbalanced data (without eliminating observations to make it balanced), specifying both variables? this way I could try to include weights = varIdent(form=~1|SoilType) instead of weights = varIdent(form=~1|ii).

Thanks.

Answer 1

There's a lot going on here so I'll take a stab at addressing some of the issues.

Unbalanced data terminology

The "unbalanced data" terminology is misleading in my view. When we have several groups and measure the value of a response variable (e.g., pH) for different numbers of study units within each group, we say that the groups are unbalanced. In your case, the groups would be all possible combinations of SoilType and LandUse. It's not clear from your description what your study unit is (e.g., soil sample) but hopefully you have several study units per group.

Combinations of levels of SoilType and LandUse

What I think you mean by "unbalanced data" is the fact that you did not consider all possible combinations of levels of SoilType and LandUse in your study - you only considered a subset of these combinations. This would explain the error message you get from R regarding the singular gls fit. You can verify if this is the case for your data using the R command:

with(mydata, table(SoilType, LandUse))

The command will show you how often each possible combination of levels of SoilType and LandUse appears in your data. If a particular combination does not appear at all in your data, the contingency table produced by the above command will show an entry of 0 for it.

Addressing missing combinations of levels of SoilType and LandUse

If you are missing some of all possible combinations of SoilType and LandUse but you try to fit a gls model including an interaction between SoilType and LandUse, you will get the dreaded singular fit message from R.

This is because the design matrix for your model is degenerate. See this excellent post on Degenerate design matrices for more details: https://rstudio-pubs-static.s3.amazonaws.com/6311_a09169ad892f4f5499874751a5fa822d.html.

To circumvent this message, you would need to "manually" define a variable which indicates what specific combination of levels of SoilType and LandUse you are dealing with for each of the study units represented in your data. For example:

mydata$Combination <- interaction(mydata$SoilType, 
                                  mydata$LandUse)

mydata$Combination <- droplevels(mydata$Combination)

Then you can fit a gls model which does not use the weights = option like this:

gls(pH ~ Combination, data = mydata) 

Here, I assume for simplicity there are no missing values in your data.

Simple example

The link on degenerate design matrices includes a simple example where Combination would be obtained by listing those combinations present in the data for two factors: f1 and f2. (In your case, f1 would be SoilType and f2 would be LandUse). In the example, f1 has 2 possible levels: a and b; f2 has 2 possible levels: A and B. All possible combinations of levels of f1 and f2 would be: aA, aB, bA, bB. However, in the data underlying that example, only the combinations aA, bA and bB are represented. If the gls() model is fitted as:

gls(y ~ 1 + Combination, data = etc.)

then R will look at the levels of Combination: aA, bA and bB. It will then set aside the first level, aA, and report 3 lines of output. The first line of output will report the estimated mean of y when Combination = aA. The second line of output will report the difference in the estimated mean value of y between Combination = bA and Combination = aA. The third line of output will report the difference in the estimated mean value of y between Combination = bB and Combination = aA. (Combination is denoted by f12 in the example).

If you compare the model gls(y ~ 1 + Combination, data = etc.) against the model gls(y ~ 1, data = etc.), you will be able to test these hypotheses:

Ho: there is no difference in the mean value of y across any of the combinations of levels of f1 and f2 represented in the study

Ha: there is a difference in the mean value of y across at least two of the combinations of levels of f1 and f2 represented in the study

The null hypothesis Ho corresponds to the model gls(y ~ 1, data = etc.), whereas the alternative hypothesis Ha corresponds to the model gls(y ~ 1 + Combination, data = etc.).

Of course, if Ho is rejected in favour of Ha, one would have to detect which specific combinations are different from each other by performing post-hoc comparisons of combinations.

Another possibility here would be that we know upfront we are only interested in comparing certain combinations - in other words, we would have some planned comparisons of combinations in mind from the study design phase: say, we would only want to compare the mean value of y for combinations bB and bA. Then we wouldn't need to test Ho against Ha - we would simply proceed to perform the comparisons of combinations of interest. This can be done by setting up and testing/estimating contrasts of combinations.

This is simple enough to do with just 2 levels per factor. But in your case, you'll have a large number of combinations, so you may either want to consolidate your LandUse definition to a smaller number of levels, or identify some specific combinations of interest, etc. Otherwise, you will have a massive number of comparisons on your hands. In any event, adjusting for multiplicity of comparisons will be a must!

Non-Constant Variability

After fitting the model gls(pH ~ Combinations, data = mydata), you can plot its residuals against SoilType to see whether there is evidence of non-constant variability across SoilType. If there is, fit a model like:

gls(pH ~ 1 + Combination, 
    weights = varIdent(form=~1|SoilType), 
    data = mydata)

Then, if you are able to fit the model without running into memory allocation issues, you would still need to consider performing multiple comparisons of combinations of levels. In particular, should you perform all possible comparisons or just a select few identified at the study planning stage?

If you can't fit the above model due to memory allocation issues, you may have no choice but to consolidate the levels of LandUse. Just use a command like:

table(mydata$LandUse)

and merge levels which are minimally represented in the data with levels that are maximally represented provided that makes substantive sense.