Deming regression finds a regression line that relates x and y variables for how they covary, i.e., it finds the relationship between x and y, whereas, for example, ordinary least squares would find the least error prediction of y given x. Does anyone have Mathematica code for Deming regression?
4
-
$\begingroup$ demonstrations.wolfram.com/… $\endgroup$– CoolwaterSep 13, 2016 at 4:02
-
$\begingroup$ @Coolwater That's nice, but perpendicular regression does not select the actual regression angles needed for Deming regression. To say that another way, if one were to rescale the y-axis variable so that the regression errors had the same x as y variances, then perpendicular regression and Deming regression would be the same. $\endgroup$– CarlOct 26, 2016 at 21:02
-
$\begingroup$ Two related questions. $\endgroup$– J. M.'s eventual burnout ♦Jan 10, 2021 at 3:09
-
$\begingroup$ @J.M.'sennui The relationship is vague for those, because they are not exactly bivariate. A more direct, but still not the same relationship would between Deming and Passing-Bablok, or Theil regressions, which unlike Deming regression are nonparametric, i.e., they use median (or ranked) slopes of all possible slopes of binary subsets of data points. $\endgroup$– CarlJan 10, 2021 at 4:06
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
No one has answered for a long time. Here is some code with Deming in it that I translated and played with from other sources. I am sorry that the annotation is not more complete. I am also not used to this site, so anyone who wants to improve the appearance or include comments in the code, go ahead.
(*Deming regression*)
n=100;
z=Table[i,{i,1,100}];
sdx=5;sdy=10;vr=(sdy/sdx)^2
alphalist={};
betalist={};
sigmaylist={};
sigmaxlist={};
nboot=10000;
alfa=0.05;
(*x=Flatten[Transpose[Table[z[[i]]+RandomVariate[NormalDistribution[0,sdx],1],{i,1,n}]]]
y=Flatten[Transpose[Table[z[[i]]+RandomVariate[NormalDistribution[0,sdy],1],{i,1,n}]]]*)
x={3.603300194803266,-0.6800782805936687,12.12247659909247,0.5720268530610206,4.043495677295204,9.548534615002781,11.60817087995439,2.976102756694152,19.931461997086238,4.05457629315898,2.3150528235849883,13.05200868281086,8.33609014173149,8.649374693340757,11.588680833263528,18.415854341788005,13.484856350672057,14.730649671993943,16.159811296238068,28.542590601303694,26.40269686158906,27.358834788118568,27.570935360098503,23.90253916277712,28.95791446667111,24.479102310685807,21.43642978054541,25.198924847189755,26.613671068480816,26.342601173804287,25.035851619716162,36.40418275093935,34.98469454643016,36.62647239823896,40.68163698254326,34.13943742666708,28.970100169689694,37.697935344366826,43.2103089112359,45.13325591572377,45.67674385080898,33.1400448483919,30.51285874093289,49.17468914927296,45.558566248123775,52.04115648231778,52.13267329178281,42.79170438181535,48.87786311460176,43.62427861082064,52.845966736399376,53.44980131769242,51.74712993517206,53.2119395788414,52.92358669132119,61.506071790666304,56.88725817597634,60.35012547108243,61.48327998754642,63.88620368163492,57.577182285642536,52.918756537899704,63.73853625711925,69.93528606880908,62.67746044369196,73.26614850129155,69.55458905854316,66.86712689530982,54.65231621751393,67.80804818705899,74.98743249205296,69.63206546301713,81.84788639291494,70.19037096082724,72.62271947368633,75.34801705387817,76.19231926597548,80.82295679328665,81.54350235028132,75.07313640705064,77.67872509988733,85.87925976532622,87.72443621581223,87.74847537101229,81.54710127589041,82.96936018691346,94.1106384321072,89.097197010559,88.92685244181976,94.25776832998227,95.99717996924906,88.45207578880714,98.41346263650759,92.96497295663569,84.44338380155591,95.19229215994184,93.89908199493311,96.94076524154501,110.17592918135131,103.49885316543727};
y={-1.8395511656654375,8.504304996609061,-19.966116834539893,25.975172366028545,21.369494172747892,-4.877412617566026,-14.172058774182155,10.145780299040053,-0.5917209901068325,12.675255681784442,-0.9869792003102731,25.040311371302707,24.606065983521926,19.820263305081582,3.630438635617173,15.63626366014385,18.394283523867287,43.206585625614835,23.107930316387687,17.839693609700003,13.619807751652411,39.47135968911134,20.94333761372782,34.12956082161715,19.82316695107436,36.53380634091501,24.830637050265786,44.69959596916163,32.341742294813095,20.92256314979675,43.03687261757748,21.849517622515023,48.80553611039267,40.79864086832398,35.91932639650507,42.23027912171956,12.687747322220925,50.16891898745286,32.734915539074805,53.9419592326429,38.52135043588691,34.297945893052024,55.68788039268125,31.512746881568745,60.062964042940706,51.66182907068255,30.846917187976242,36.179390424392274,49.446482305660034,50.49434199541827,62.04749084328972,26.386935795317644,58.9981128144914,67.92048978928926,62.7506803556757,62.39259624435499,66.5434886534815,67.78924742545377,62.82849811036929,73.64821434484001,64.72200573642179,60.96313872454566,69.32972260256984,53.77917194358318,66.50672596000553,76.18725716807197,60.03238207075945,69.90527731072103,64.23587954374602,63.49569676235791,71.2145386216645,91.18862163690152,61.26540465241208,55.719076445687435,76.24217744645262,77.0495478763569,69.98869923448973,73.979994284022,86.75489711361934,69.32326936982429,67.73890110672403,84.69847299952067,85.1083933762478,84.15598167973721,102.57939622396309,69.94947069853325,87.25047997362489,86.5225455632367,100.73914245483502,91.23708149430198,81.60545270046555,100.09422669364191,104.55087400322526,73.61619077924142,87.21761388748021,106.54065896031568,103.59361698169805,98.914445398764,111.38420843501999,92.53700682734804};
SSDy=Variance[y]*(n-1)
SSDx=Variance[x]*(n-1)
SPDxy=Covariance[x,y]*(n-1)
beta0=beta=(SSDy-vr*SSDx+Sqrt[(SSDy-vr*SSDx)^2+4*vr*SPDxy^2])/(2*SPDxy);
alpha0=alpha=Mean[y]-Mean[x]*beta;
Print["y = ",alpha," + ",beta," x"]
ksi[i_]:=(vr*x[[i]]+beta*(y[[i]]-alpha))/(vr+beta^2)
sigmax=N[1/((n-2)*vr) (vr*Sum[(x[[i]]-ksi[i])^2,{i,1,n}]+Sum[(y[[i]]-alpha-beta*ksi[i])^2,{i,1,n}])]
sigmay=vr*sigmax
sigmax=Sqrt[sigmax]
sigmay=Sqrt[sigmay]
(*mx=Mean[x]
my=Mean[y]
Print["SDm = ",SDm=Sqrt[Sum[(y[[i]]-my)^2,{i,1,n}]/(n-2)]/Sqrt[Sum[(x[[i]]-mx)^2,{i,1,n}]]]
Needs["HypothesisTesting`"]
StudentTCI[beta, SDm/n^(1/2)/vr,n-1]
StudentTCI[beta, 4 SDm/n^(1/2),n-1]
InverseCDF[StudentTDistribution[beta, 4 SDm/n^(1/2),n-1],alfa/2]*)
plot=Plot[alpha+beta z,{z,Min[x],Max[x]}];
lplot=ListPlot[Table[{x[[i]],y[[i]]},{i,1,n}]];
dat=Table[{x[[i]],y[[i]]},{i,1,n}];
ols=Fit[dat,{1,w},w]
plotOLS=Plot[ols,{w,Min[x],Max[x]},PlotStyle->Red];
Show[lplot,plot,plotOLS]
(*alphalist={alpha};
betalist={beta};
sigmaylist={sigmay};
sigmaxlist={sigmax};*)
(*bootstrap*)
Do[new=RandomChoice[dat,n];
x=new[[All,1]];
y=new[[All,2]];
SSDy=Variance[y]*(n-1);
SSDx=Variance[x]*(n-1);
SPDxy=Covariance[x,y]*(n-1);
beta=(SSDy-vr*SSDx+Sqrt[(SSDy-vr*SSDx)^2+4*vr*SPDxy^2])/(2*SPDxy);
alpha=Mean[y]-Mean[x]*beta;
(*Print["y = ",alpha," + ",beta," x"];*)
ksi[i_]:=(vr*x[[i]]+beta*(y[[i]]-alpha))/(vr+beta^2);
sigmax=N[1/((n-2)*vr) (vr*Sum[(x[[i]]-ksi[i])^2,{i,1,n}]+Sum[(y[[i]]-alpha-beta*ksi[i])^2,{i,1,n}])];
sigmay=vr*sigmax;
sigmax=Sqrt[sigmax];
sigmay=Sqrt[sigmay];
alphalist=AppendTo[alphalist,alpha];
betalist=AppendTo[betalist,beta];
sigmaylist=AppendTo[sigmaylist,sigmay];
sigmaxlist=AppendTo[sigmaxlist,sigmax],{count,1,nboot}]
Print["Mean, SD of \[Alpha] = ",m\[Alpha]=Mean[alphalist],", ",\[Sigma]\[Alpha]=StandardDeviation[alphalist],"\n","Mean, SD of \[Beta] = ",m\[Beta]=Mean[betalist],", ",\[Sigma]\[Beta]=StandardDeviation[betalist],"\n",
"Mean, SD of Subscript[\[Sigma], y] = ",Mean[sigmaylist],", ", StandardDeviation[sigmaylist],"\n",
"Mean, SD of Subscript[\[Sigma], x] = ",Mean[sigmaxlist],", ", StandardDeviation[sigmaxlist]]
Print[Histogram[alpha0-alphalist,ImageSize->Small],Histogram[betalist-beta0+1,ImageSize->Small]]
Print["95% CI \[Alpha] ",Qalphalow=Quantile[alphalist-alpha0,0.025,{{0,1},{0,1}}]," to ",
Qalphahigh=Quantile[alphalist-alpha0,0.975,{{0,1},{0,1}}],"\tmean = ",1/2Qalphalow+1/2Qalphahigh,"\tmedian = ",Median[alphalist-alpha0],"\trange = ",Qalphahigh-Qalphalow]
Print["95% CI \[Beta] ",Qbetalow=Quantile[betalist-beta0+1,0.025,{{0,1},{0,1}}]," to ",
Qbetahigh=Quantile[betalist-beta0+1,0.975,{{0,1},{0,1}}],"\tmean = ",1/2Qbetalow+1/2Qbetahigh,"\tmedian = ",Median[betalist-beta0+1]"\trange = ",Qbetahigh-Qbetalow]
Print["emperical p(\[Alpha]) = 0 is ",empalpha=NumberForm[1-2 Abs[1/2-CDF[EmpiricalDistribution[alphalist],0]],Round[Log10[nboot]]]]
Print["emperical p(\[Beta]) = 1 is ",empbeta=NumberForm[1-2 Abs[1/2-CDF[EmpiricalDistribution[betalist],1]],Round[Log10[nboot]]]]
OUTPUT
4
92078.
87071.2
82188.5
y = 2.02335 + 0.976973 x
30.0562
120.225
5.48236
10.9647
3.69673 +0.943923 w
Mean, SD of \[Alpha] = 1.98257, 2.90811
Mean, SD of \[Beta] = 0.977749, 0.0447706
Mean, SD of Subscript[\[Sigma], y] = 10.8317, 0.749437
Mean, SD of Subscript[\[Sigma], x] = 5.41583, 0.374718
95% CI \[Alpha] -5.91613 to 5.56346 mean = -0.176334 median = 0.0300523 range = 11.4796
95% CI \[Beta] 0.914103 to 1.0917 mean = 1.0029 median = 1.00034 range = 0.177593
emperical p(\[Alpha]) = 0 is 0.486
emperical p(\[Beta]) = 1 is 0.6034
$\begingroup$
$\endgroup$
3
See "Least Squares Fitting-Perpendicular Offsets" by Dimitris Apostolos Sardelis & Theodoros Valahas
DOI: 10.13140/RG.2.2.24224.84488 (available from Research Gate - see (3) below)
(1) The current article first appeared at the Wolfram Library Archive: Dimitris Sardelis and Theodore Valahas: “Least SquaresFitting-Perpendicular Offsets.” (2004) http://library.wolfram.com/infocenter/MathSource/5292/ Its present form is mainly a re-edition effected by the use of Mathematica, version 11.2.
(2) Eric W. Weisstein: "Least Squares Fitting--Perpendicular Offsets." (2002) From Wolfram MathWorld--A Wolfram WebResource. http://mathworld.wolfram.com/LeastSquaresFittingPerpendicularOffsets.html6 Least Squares Fitting-Perpendicular Offsets (Re-edited).nb
(3) (PDF) Least Squares Fitting-Perpendicular Offsets. Available from: https://www.researchgate.net/publication/315069657_Least_Squares_Fitting-Perpendicular_Offsets [accessed Oct 20 2023].
A Mathematica Example is included
-
2$\begingroup$ Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. $\endgroup$– Community BotOct 19 at 0:17
-
-
$\begingroup$ Perpendicular least squares is not the same as Deming regression. In specific, the regression angle for Deming regression is a variable, not a right angle. $\endgroup$– CarlOct 20 at 3:16