Drift correction for geochemical data

Question

Emily 2023-6-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1984369-drift-correction-for-geochemical-data

评论： Star Strider 2023-6-16

fiji.csv

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-6-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1984369-drift-correction-for-geochemical-data#answer_1257634

在 MATLAB Online 中打开

fiji.csv

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.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2023-6-16

在 MATLAB Online 中打开

fiji.csv

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 2023-6-16

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 2023-6-16

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.)

请先登录，再进行评论。