12
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ demonstrations.wolfram.com/… $\endgroup$
    – Coolwater
    Sep 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$
    – Carl
    Oct 26, 2016 at 21:02
  • $\begingroup$ Two related questions. $\endgroup$ 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$
    – Carl
    Jan 10, 2021 at 4:06

2 Answers 2

8
$\begingroup$

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

enter image description here

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

enter image description hereenter image description here

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
$\endgroup$
2
$\begingroup$

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

$\endgroup$
3
  • 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 Bot
    Oct 19 at 0:17
  • $\begingroup$ Provided all the information available $\endgroup$ Oct 20 at 1:22
  • $\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$
    – Carl
    Oct 20 at 3:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.