10
$\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$
  • $\begingroup$ demonstrations.wolfram.com/… $\endgroup$ – Coolwater Sep 13 '16 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 '16 at 21:02
4
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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