Drift correction for geochemical data

I'm working on a code to automate data processing. Towards the end of my code, I need to determine if isotope data has drifted throughout the course of the analysis. Data are attached.

The first column of the data – fiji(:,1) – shows the line number of the analysis (each measurement is a new line)

The fourth column of the data – fiji(:,4) – shows the geochemical value for each measurement

What I need help with is a succinct way to determine if the slope of the drift is significantly different from zero. If the slope is not different from zero, then we don't need to worry about drift. If it is different from zero, I'll look at the R^2 value to determine if I want to proceed from a drift correction (already included in the code_.

Any insight about coding up if the slope is different from zero would be very appreciated. My first instinct was to use ANOVA, but I'm not sure how to apply it to this dataset.

Thanks!!

fiji = readtable("fiji.csv");
fiji = table2array(fiji);

%Calculate drift correction
drift_d18O = fitlm(fiji(:,1), fiji(:,4));
drift_d18O_slope = drift_d18O.Coefficients{2,1};
drift_d18O_intercept = drift_d18O.Coefficients{1,1};

%%If drift slope is significantly different from zero, then check R^2 below
%%If drift not significantly different from zero, stop analysis here & skip
%%to end

if (drift_d18O.Rsquared.Adjusted > 0.7)
    disp('R squared indicates there IS significant d18O drift')
else 
    disp('R squared indicates there IS NOT significant d18O drift')
end 

%Plot the drift correction
figure(1)
scatter(fiji(:,1), fiji(:,4));

%drift functions

function d = drift_d18O_func(x, d18O_slope, d18O_intercept)
    d = x.*d18O_slope+d18O_intercept;
end

Answer 1

[Star Strider]

First, plot the data —

fiji = readmatrix('fiji.csv')

fiji = 12×8

1.0e+04 * 0.0040 0.0007 0.0004 -0.0005 -0.0026 1.6368 0.0000 0.0001 0.0041 0.0007 0.0005 -0.0005 -0.0026 1.6379 0.0000 0.0001 0.0042 0.0007 0.0006 -0.0005 -0.0026 1.6377 0.0000 0.0001 0.0178 0.0007 0.0004 -0.0005 -0.0025 1.6268 0.0000 0.0001 0.0179 0.0007 0.0005 -0.0005 -0.0025 1.6169 0.0000 0.0001 0.0180 0.0007 0.0006 -0.0005 -0.0025 1.6244 0.0000 0.0001 0.0316 0.0007 0.0004 -0.0005 -0.0025 1.6219 0.0000 0.0001 0.0317 0.0007 0.0005 -0.0005 -0.0025 1.6212 0.0000 0.0001 0.0318 0.0007 0.0006 -0.0005 -0.0025 1.6093 0.0000 0.0001 0.0454 0.0007 0.0004 -0.0005 -0.0024 1.6252 0.0000 0.0001

Col1 = fiji(:,1)

Col1 = 12×1

40 41 42 178 179 180 316 317 318 454

Col4 = fiji(:,4)

Col4 = 12×1

-5.0754 -5.0321 -5.0472 -5.0311 -5.0261 -5.1097 -4.9596 -5.1198 -5.0885 -5.1389

% d4d1 = gradient(fiji(:,4)) ./ gradient(fiji(:,1));

figure
plot(fiji(:,1), fiji(:,4), '.-')
xlabel('Col #1')
ylabel('Col #4')
grid

What do you want to get from these data?

I doubt that a linear fit will provide any useful information.

.

[Emily]

Some clarification –

This code will be used to process many data files. Sometimes there is a very strong linear trend where Col 4 changes over Col 1 (in the case of this file, there doesn't appear to be a very strong correlation btwn Col 4 and Col 1).

But, I need to figure out how to code a succinct way to figure out (for each data file) if Col 4 is changing with Col 1. If it does change, then I need to apply a drift correction. I use a linear fit btwn Col 1 and Col 4 to accomplish this. If Col 4 does not change significantly with Col 1, I do not need to apply a drift correction.

What I would like to do is be able to see if the slope of the linear fit between Col 4 and Col 1 is significantly different from zero, as I think this will be the best way to tell me if I need to do a drift correction or not. I think some kind of test like anova might be the way to go, but this is what I would like insight on!

[Star Strider]

I would use the slope probability instead, since it indicates exactly what you want to determine —

fiji = readtable("fiji.csv");
fiji = table2array(fiji);

%Calculate drift correction
drift_d18O = fitlm(fiji(:,1), fiji(:,4))

drift_d18O = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                    Estimate          SE         tStat       pValue  
                   ___________    __________    _______    __________

    (Intercept)        -5.0205      0.032756    -153.27    3.4335e-18
    x1             -0.00024866    0.00011215    -2.2172      0.050937


Number of observations: 12, Error degrees of freedom: 10
Root Mean Squared Error: 0.0599
R-squared: 0.33,  Adjusted R-Squared: 0.263
F-statistic vs. constant model: 4.92, p-value = 0.0509

drift_d18O_slope = drift_d18O.Coefficients{2,1};
drift_d18O_intercept = drift_d18O.Coefficients{1,1};

%%If drift slope is significantly different from zero, then check R^2 below
%%If drift not significantly different from zero, stop analysis here & skip
%%to end

if (drift_d18O.Coefficients.pValue(2) > 0.05)
    disp('Slope p-value indicates there IS significant d18O drift')
else 
    disp('Slope p-value indicates there IS NOT significant d18O drift')
end 

Slope p-value indicates there IS significant d18O drift

xv = linspace(min(fiji(:,1)), max(fiji(:,1))).';
[y yci] = predict(drift_d18O, xv);

%Plot the drift correction
figure(1)
hs = scatter(fiji(:,1), fiji(:,4), 'DisplayName','Data');
hold on
hr = plot(xv, y, '-r', 'DisplayName','Regression');
hc = plot(xv, yci, '--r', 'DisplayName','95% Confidence Intervals');
hold off
legend([hs hr hc(1)], 'Location','best')

%drift functions

function d = drift_d18O_func(x, d18O_slope, d18O_intercept)
    d = x.*d18O_slope+d18O_intercept;
end

.

[Emily]

Thank you! I may be misunderstanding but I think it should be the reverse? If the p value < 0.05 then the drift IS significant because the slope is significant?

[Star Strider]

My pleasure!

If the p-statistic is less than 0.05, the slope is significantly different from zero, so the inequality test should be reversed in that test. (I didn’t read the code carefully enough.)