# Simple regression on subsets vs ANCOVA with interaction term

I'm attempting to predict a continuous variable as a function of one continuous variable and one categorical variable. This led me to an ANCOVA, but I've violated the homogeneity of regression slopes assumption. My next thought was to simply subset my data by the categorical variable and fit a simple regression to each subset, but this feels problematic as well. When I attempted this in R, I received identical estimates for each term but wildly different SEs (and therefore significance values). I thought that these two methods would produce the same output because the model matrix is essentially performing subsetting with its zeroes and ones, but now I'm wondering which of the two is more accurate to report.

My main concern is assessing the confidence of predictions made with this model. If someone gave me these two predictors and asked for the predicted output, the models would agree - but if they asked for a confidence estimate I'd have no idea which of the two outputs to report. I assume the ANCOVA output is questionable because of the assumption violation, but the perfect predictor agreement makes me doubt whether it's actually unsuitable.

Reproducible example in R below:

library(dplyr)
nit_data <- structure(
list(no23_um = c(
0.0015, 0.7835, 0.3855, 0.0035, 0.2708, 2.2198, 0.0032, 0.4371, 6.0165,
0.0032, 0.0187, 0.0036, 0.002, 0.8623, 2.6331, 1.0191, 5.2738, 2.4266,
0.0021, 0.007, 0.4422, 0.0031, 0.1195, 2.7913, 0.0026, 4.4071, 0.6548,
0.0043, 0.1485, 0.0041, 0.0019, 0.8997, 0.1487),
sla = c(
-0.0392, 0.1837, 0.252, 0.228, 0.1793, 0.1036, 0.0111, -0.0883,
-0.1609, -0.0685, 0.073, -0.0392, -0.0883, -0.1609, -0.0692, 0.073,
-0.0392, 0.1837, 0.1837, 0.252, 0.228, 0.1793, 0.1036, 0.0111,
0.1036, -0.0883, -0.0692, 0.252, 0.228, 0.073, -0.1609, 0.1793, 0.0111),
abs_depth = c(
15, 122, 175, 15, 122, 175, 15, 110, 175, 15, 105, 100, 15, 103,
175, 175, 175, 175, 15, 119, 175, 15, 110, 175, 15, 175, 112,
15, 133, 15, 15, 175, 110)
), row.names = c(NA, -33L), class = "data.frame") %>%
mutate(depth=case_when(
abs_depth==15~"Surface",
abs_depth==175~"Deep",
TRUE~"DCM"
)) %>%
mutate(depth=factor(depth, levels = c("Surface", "DCM", "Deep")))

# ANCOVA with interaction:
broom::tidy(lm(no23_um~depth*sla, data=nit_data))

# Subset and simple regression:
lapply(c("Surface", "DCM", "Deep"), function(x){
broom::tidy(lm(no23_um~sla, data=nit_data[nit_data$depth==x,])) %>% mutate(term=paste0(x, ":", term)) }) %>% bind_rows() %>% slice(c(1,3,5,2,4,6))  Output: > # ANCOVA with interaction: > broom::tidy(lm(no23_um~depth*sla, data=nit_data)) # A tibble: 6 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 0.00263 0.201 0.0131 9.90e- 1 2 depthDCM 0.371 0.285 1.30 2.04e- 1 3 depthDeep 3.32 0.285 11.7 4.74e-12 4 sla 0.00376 1.37 0.00274 9.98e- 1 5 depthDCM:sla -0.975 1.94 -0.503 6.19e- 1 6 depthDeep:sla -12.0 1.94 -6.18 1.33e- 6 > # Subset and simple regression: > lapply(c("Surface", "DCM", "Deep"), function(x){ + broom::tidy(lm(no23_um~sla, data=nit_data[nit_data$depth==x,] .... [TRUNCATED]
# A tibble: 6 x 5
term                 estimate std.error statistic    p.value
<chr>                   <dbl>     <dbl>     <dbl>      <dbl>
1 Surface:(Intercept)   0.00263  0.000264      9.98 0.00000363
2 DCM:(Intercept)       0.373    0.101         3.69 0.00502
3 Deep:(Intercept)      3.32     0.334         9.96 0.00000369
4 Surface:sla           0.00376  0.00180       2.09 0.0658
5 DCM:sla              -0.971    0.689        -1.41 0.192
6 Deep:sla            -12.0      2.27         -5.27 0.000514