Hoofdstuk 9 Measurement Error Correction

Het R package mecor (Nab et al., 2021; Nab, 2021) bevat de mogelijkheid om de data uit een substudie te gebruiken om te corrigeren voor measurement error in de predictor (of in de outcome, maar dat is voor het SoilHarmony project nvt). De functies uit dit package geven zowel de niet-gecorrigeerde (proxy relatie) als de gecorrigeerde regressie (functionele relatie). De niet-gecorrigeerde regressie kan gebruikt worden voor predicties te maken (transferfunctie) en de gecorrigeerde regressie schat de functionele relatie (beide interessant).

# ter vergelijking met de mecor-methoden printen we nog eens de relatie 
# Y vs. X:
(lm(pHH2O ~ pHCa, data = data_wosis_pH))

## 
## Call:
## lm(formula = pHH2O ~ pHCa, data = data_wosis_pH)
## 
## Coefficients:
## (Intercept)         pHCa  
##      1.0425       0.9118

# Y vs. W:
(lm(pHH2O ~ pHCa_obs, data = data_wosis_pH))

## 
## Call:
## lm(formula = pHH2O ~ pHCa_obs, data = data_wosis_pH)
## 
## Coefficients:
## (Intercept)     pHCa_obs  
##      1.7455       0.7957

9.1 Replicatie-substudie

De substudie met herhaalde metingen van de predictor kan als volgt gebruikt worden in de mecor functie. We gebruiken maximum-likelihood estimation voor het model te fitten (bij de standaard methode zijn de resultaten afhankelijk van de volgorde waarin de replicaten staan). De schatting van de gecorrigeerde regressie-coëfficiënt ligt zeer dicht bij de waarde uit de volledige dataset (zie hierboven) en het 95% CI omvat de waarde. Deze methode lijkt dus beter dan de functies uit het mcr en smatr package, alvast voor deze dataset.

mecor_m1 <- mecor(
  pHH2O ~ MeasError(pHCa_obs, replicate = cbind(pHCa_obs2, pHCa_obs3)),
  method = "mle", 
  data = data_train
)
summary(mecor_m1)

## 
## Call:
## mecor(formula = pHH2O ~ MeasError(pHCa_obs, replicate = cbind(pHCa_obs2, 
##     pHCa_obs3)), data = data_train, method = "mle")
## 
## Coefficients Corrected Model:
##              Estimate       SE
## (Intercept)  1.055635 0.044965
## cor_pHCa_obs 0.909343 0.007405
## 
## 95% Confidence Intervals:
##              Estimate      LCI      UCI
## (Intercept)  1.055635 0.967504 1.143766
## cor_pHCa_obs 0.909343 0.894829 0.923856
## 
## The measurement error is corrected for by application of maximum likelihood estimation 
## 
## Coefficients Uncorrected Model:
##              Estimate Std. Error t value  Pr(>|t|)
## (Intercept) 1.7430264  0.0117642  148.16 < 2.2e-16
## pHCa_obs    0.7960465  0.0018823  422.90 < 2.2e-16
## 
## 95% Confidence Intervals:
##             Estimate      LCI      UCI
## (Intercept) 1.743026 1.719968 1.766085
## pHCa_obs    0.796047 0.792357 0.799736
## 
## Residual standard error: 0.5393357 on 34871 degrees of freedom

Het niet-gecoorigeerde regressie-model kan gebruikt worden voor predicties te maken. Hieronder staan de predicties met coverage en RMSE.

Figure 9.1: Evaluatie van de test-dataset met de transferfunctie.

Extra covariaten kunnen gebruikt worden in het mecor package (er wordt veronderstelt dat deze zonder meetfout zijn gemeten). De RMSE wordt hierdoor kleiner.

mecor_m2 <- mecor(
  pHH2O ~ MeasError(pHCa_obs, replicate = cbind(pHCa_obs2, pHCa_obs3)) +
    bulk_dens,
  method = "mle", 
  data = data_train
)
summary(mecor_m2)

## 
## Call:
## mecor(formula = pHH2O ~ MeasError(pHCa_obs, replicate = cbind(pHCa_obs2, 
##     pHCa_obs3)) + bulk_dens, data = data_train, method = "mle")
## 
## Coefficients Corrected Model:
##              Estimate       SE
## (Intercept)  0.923218 0.039330
## cor_pHCa_obs 0.906203 0.007580
## bulk_dens    0.098382 0.012499
## 
## 95% Confidence Intervals:
##              Estimate      LCI      UCI
## (Intercept)  0.923218 0.846132 1.000304
## cor_pHCa_obs 0.906203 0.891347 0.921060
## bulk_dens    0.098382 0.073885 0.122879
## 
## The measurement error is corrected for by application of maximum likelihood estimation 
## 
## Coefficients Uncorrected Model:
##              Estimate Std. Error t value  Pr(>|t|)
## (Intercept) 1.4669134  0.0187639  78.177 < 2.2e-16
## pHCa_obs    0.7904951  0.0018959 416.941 < 2.2e-16
## bulk_dens   0.2012385  0.0106886  18.827 < 2.2e-16
## 
## 95% Confidence Intervals:
##             Estimate      LCI      UCI
## (Intercept) 1.466913 1.430135 1.503691
## pHCa_obs    0.790495 0.786779 0.794211
## bulk_dens   0.201238 0.180288 0.222188
## 
## Residual standard error: 0.5366228 on 34870 degrees of freedom

Figure 9.2: Evaluatie van de test-dataset met de transferfunctie.

9.2 Validatie-substudie (intern)

De substudie waarbij de predictor zonder meetfout wordt gemeten kan als volgt gebruikt worden in de mecor functie. Ook hier kunnen extra covariaten gebruikt worden. We gebruiken efficiënt estimation om het model te fitten. De schatting van de gecorrigeerde regressie-coëfficiënt ligt zeer dicht bij de waarde uit de volledige dataset (zie hierboven) en het 95% CI omvat de waarde. Deze methode lijkt dus beter dan de functies uit het mcr en smatr package, alvast voor deze dataset. De predicties van dit model hebben dezelfde coverage en RMSE als het vorige want de predicties gebeuren obv het niet-gerorrigeerde model.

mecor_m3 <- mecor(
  pHH2O ~ MeasError(pHCa_obs, reference = pHCa_sub) +
    bulk_dens,
  method = "efficient", 
  data = data_train
)
summary(mecor_m3)

## 
## Call:
## mecor(formula = pHH2O ~ MeasError(pHCa_obs, reference = pHCa_sub) + 
##     bulk_dens, data = data_train, method = "efficient")
## 
## Coefficients Corrected Model:
##             Estimate       SE
## (Intercept) 0.913947 0.026590
## pHCa_sub    0.910684 0.002915
## bulk_dens   0.089958 0.014766
## 
## 95% Confidence Intervals:
##             Estimate      LCI      UCI
## (Intercept) 0.913947 0.861832 0.966063
## pHCa_sub    0.910684 0.904971 0.916397
## bulk_dens   0.089958 0.061018 0.118898
## 
## The measurement error is corrected for by application of efficient regression calibration 
## 
## Coefficients Uncorrected Model:
##              Estimate Std. Error t value  Pr(>|t|)
## (Intercept) 1.4669134  0.0187639  78.177 < 2.2e-16
## pHCa_obs    0.7904951  0.0018959 416.941 < 2.2e-16
## bulk_dens   0.2012385  0.0106886  18.827 < 2.2e-16
## 
## 95% Confidence Intervals:
##             Estimate      LCI      UCI
## (Intercept) 1.466913 1.430135 1.503691
## pHCa_obs    0.790495 0.786779 0.794211
## bulk_dens   0.201238 0.180288 0.222188
## 
## Residual standard error: 0.5366228 on 34870 degrees of freedom

Figure 9.3: Evaluatie van de test-dataset met de transferfunctie.

9.3 Validatie-substudie (extern)

De validatie-substudie hierboven was intern, in de zin dat de locaties uit de substudie ook gemeten zijn in de hoofdstudie. Indien dit niet het geval is en de locaties uit een andere studie komen (maar wel tot dezelfde populatie behoren natuurlijk), kan die info ook gebruikt worden in het mecor package. We fitten daarom eerst een lineair model in de externe validatie-dataset en gebruiken dit model om het model in de hoofdstudie te calibreren. Uit de resultaten hieronder blijkt dat de correctie niet zo goed is als met de andere methoden hierboven (de gecorrigeerde \(\beta\) is verder van de nominale waarde dan voorheen), maar de predicties zijn natuurlijk nog steeds bruikbaar want die gebeuren obv het niet-gecorrigeerde model.

# gebruik de validatie-substudie als was het een externe studie
ext_lm <- lm(
  pHH2O ~ pHCa_sub + bulk_dens,
  data = data_train %>% filter(substudy == 1)
)

mecor_m4 <- mecor(
  pHH2O ~ MeasErrorExt(pHCa_obs, model = ext_lm) + bulk_dens,
  data = data_train
)
summary(mecor_m4)

## 
## Call:
## mecor(formula = pHH2O ~ MeasErrorExt(pHCa_obs, model = ext_lm) + 
##     bulk_dens, data = data_train)
## 
## Coefficients Corrected Model:
##              Estimate       SE
## (Intercept)  0.672841 0.034111
## cor_pHCa_obs 0.867696 0.003754
## bulk_dens    0.124875 0.018173
## 
## 95% Confidence Intervals:
##              Estimate      LCI      UCI
## (Intercept)  0.672841 0.605985 0.739696
## cor_pHCa_obs 0.867696 0.860339 0.875053
## bulk_dens    0.124875 0.089257 0.160493
## 
## The measurement error is corrected for by application of regression calibration 
## 
## Coefficients Uncorrected Model:
##              Estimate Std. Error t value  Pr(>|t|)
## (Intercept) 1.4669134  0.0187639  78.177 < 2.2e-16
## pHCa_obs    0.7904951  0.0018959 416.941 < 2.2e-16
## bulk_dens   0.2012385  0.0106886  18.827 < 2.2e-16
## 
## 95% Confidence Intervals:
##             Estimate      LCI      UCI
## (Intercept) 1.466913 1.430135 1.503691
## pHCa_obs    0.790495 0.786779 0.794211
## bulk_dens   0.201238 0.180288 0.222188
## 
## Residual standard error: 0.5366228 on 34870 degrees of freedom

Figure 9.4: Evaluatie van de test-dataset met de transferfunctie.

9.4 Andere substudies

Het mecor package kan ook gebruik maken van calibratie-substudies. Die heb ik nog niet geprobeerd in deze dataset. Te bekijken of dat van toepassing is voor SoilHarmony of niet.

Referenties

Nab L. (2021). Mecor: Measurement error correction in linear models with a continuous outcome. https://github.com/LindaNab/mecor.

Nab L., Van Smeden M., Keogh R.H. & Groenwold R.H.H. (2021). Mecor: An R package for measurement error correction in linear regression models with a continuous outcome. Computer Methods and Programs in Biomedicine 208: 106238. https://doi.org/10.1016/j.cmpb.2021.106238.