2

I have fit a 3-knot restricted cubic spline model to my data using the rms package in R with the following code:

newmod <- ols(Carbon ~ rcs(Transect, 3) * Location, data = datc2, x = TRUE, y = TRUE)

In my model, Carbon is the response variable, transect is a continuous explanatory variable, and location is a categorical explanatory variable with 2 levels (AL and FL). The model with Transect, Location, and a Transect*Location interaction was the best model according to AIC & BIC. When I plot the model predictions, I get the following graph:

enter image description here

The difference between the 2 areas is clear, but now what I want to figure out is at what levels of Transect are the Carbon values for areas AL and FL significantly different. Visually, the confidence intervals overlap from about transect 4 to transect 8, so it seems like there is no significant difference between the 2 Locations for those Transect values. However, I know just looking at overlap in the confidence intervals is not a hypothesis test for differences between the Locations. What is the best approach for examining the difference in Carbon between the two Locations at a given value of Transect? Would I need to perhaps bootstrap the difference between the two Locations at each Transect value? If so, how is that done in R?

FishMasterB
  • 45
  • 4

1 Answers1

2

If you want to examine the differences between the two values of Location at different Transect levels, you can use the regression coefficients returned by the model directly (provided that the model is adequately validated and calibrated).

Each estimated outcome difference from a linear model is just a linear combination of the regression coefficients, weighted by the corresponding difference in predictor values. The standard error of the difference is given by the formula for the variance of a sum of correlated variables, taking into account the (co)variances among the coefficient estimates. If you are making multiple such comparisons, you use appropriate multiple-testing correction.

As you are already using the rms package, implementation could be via its contrast.rms() function. You can specify any combination of different predictor values you wish and get the estimated outcome difference, confidence limits, etc. Play a bit with some simple cases to make sure that you have the syntax under control, as the flexibility that function provides necessarily comes with some complexity in how to use it. If you wish, you can also use that function to use bootstrapped estimates instead of the estimates based directly on the regression-coefficient estimates.

If those are 95% confidence limits, then you will have a significant Location difference for most of your Transect values except near where the curves cross. Non-overlap of 95% confidence intervals can be on the order of p < 0.005 for the corresponding difference

EdM
  • 57,766
  • 7
  • 66
  • 187
  • Thank you, the contrast.rms function works great! Using this function, I can get a table of contrasts for each comparison, standard errors, confidence limits, t statistics, and P-values. Is it safe to say that if the t statistic is significant (p<0.05), there is a significant difference in carbon between the two locations at that transect value? – FishMasterB Jun 23 '21 at 13:27
  • @FishMasterB that's the usual criterion for statistical significance. Depending on how you did this, you might need to do a multiple-comparison correction. I think that the significant interaction is the most striking result here; just which transects happen to show a "significant" difference is more a matter of how many samples you analyzed. – EdM Jun 23 '21 at 13:32