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 lack of A.I. ♦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
|
1 Answer
$\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