runit.mcPBequi.R
In mcr: Method Comparison Regression

cat("\n\nmcPBequi.R method comparison test cases\n\n")

## check values for equiariant Passing-Bablok and Theil-Sen regression against regression coefficients computed with the deming package (v. 1.4) by Terry Therneau. 
## Confidence intervals are calculated by another method and are therefore not directly comparable. 

#library(deming)

test.mc.PBequi.call1 <- function() 
{
    data(creatinine)
    crea <- creatinine[complete.cases(creatinine),]

# Compare to deming package
#    library(deming)

#    Compare equivariant PaBa regression

#    pbdeming <- pbreg(serum.crea~plasma.crea,crea, method=3)
#    dput(pbdeming[[1]])
#Output: c(`(Intercept)` = 0.102307692307692, plasma.crea = 0.923076923076924 

    pbeqsmall <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="PBequi",calcCI=TRUE,methodlarge=F)
    pbeqlarge <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="PBequi",calcCI=TRUE,methodlarge=T)
    
    
    
    # values are copied and pasted from deming package
    
                   
    checkEquals( 0.102307692307692, pbeqsmall$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.923076923076924, pbeqsmall$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 14 decimal places
     checkEquals( 0.102307692307692, pbeqlarge$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 14 decimal places
    checkEquals( 0.923076923076924, pbeqlarge$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 14 decimal places

# Compare Theil-Sen regression

#     tsdeming <- theilsen(serum.crea~plasma.crea,crea )  
#     dput(tsdeming[[1]])   
     tssmall <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="TS",calcCI=TRUE,methodlarge=F)
    tslarge <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="TS",calcCI=TRUE,methodlarge=T)

    checkEquals( 0.17796875, tssmall$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.84375, tssmall$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
     checkEquals( 0.17796875, tslarge$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.84375, tslarge$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
  
    X.td9 <- c( 2.47945292480324E-5, 5.4492735891554E-5, 7.1400324331745E-5, 8.24418028768875E-6, 5.25094938140207E-5, 4.96071630752203E-5, 8.38275749515994E-5, 
                4.79953883994996E-5, 6.53945975719383E-5, 6.06219105809957E-5, 5.48301410574115E-6, 8.46612595956848E-5, 7.38714176242463E-5, 4.79437559865663E-5, 
                5.96173651220441E-5, 5.97517181536454E-5, 6.75262023095708E-5, 4.37996746882773E-5, 4.5633493334196E-5, 2.67327933066674E-5, 5.39428135962028E-5, 
                7.80769623282439E-5, 9.68983239438709E-5, 7.70737197541118E-5, 2.91866029770011E-5, 3.95302214111076E-5, 7.82300704374343E-5, 6.8911992683726E-5, 
                1.35164373215798E-5, 2.82443848442956E-5, 4.67350708937759E-5, 7.97870204301856E-5, 2.61963488304632E-5, 6.72743001782288E-5, 5.90586328684933E-5, 
                3.43475612827838E-5, 7.07269059209775E-5, 1.38876686310818E-5, 1.52437947219811E-5, 5.17623575861341E-5, 3.77753027538479E-6, 4.45374962805357E-5, 
                7.32847478866719E-5, 2.41852586285496E-5, 1.44586568718652E-5, 3.15521957485487E-5, 4.04157159760381E-6, 7.24707888984188E-5, 2.88848969586527E-5, 
                1.81454715994012E-5 )
    
    Y.td9 <- c( 2.11052371498578E-5, 4.23169020455709E-5, 6.52087072020304E-5, 9.24793174778137E-6, 6.66343609634913E-5, 3.79576987793361E-5, 9.56084903375025E-5, 
                5.78813367248612E-5, 5.9825287298063E-5, 6.29197167347258E-5, 6.45317326688806E-6, 6.81545583186772E-5, 6.14328864897738E-5, 3.62479512864404E-5, 
                4.53666706876532E-5, 4.51337596066822E-5, 5.88801541326545E-5, 5.35866863682414E-5, 5.10872483563812E-5, 2.00470787909169E-5, 5.52967198327848E-5, 
                7.34531006344544E-5, 8.56661460292695E-5, 6.9278207448691E-5, 2.77691478340553E-5, 4.50792600802043E-5, 7.4091530765788E-5, 8.51728132968523E-5, 
                1.48162660519741E-5, 2.81423722618154E-5, 4.70798411549798E-5, 9.26814606305503E-5, 2.7362811205277E-5, 6.75935887363532E-5, 5.91235046402203E-5, 
                4.20744070566868E-5, 5.82538662865107E-5, 1.15176888782602E-5, 1.37798955448975E-5, 5.95000570904813E-5, 3.87235051249896E-6, 3.59395060007184E-5, 
                0.000100595978033436, 1.93943506560089E-5, 1.27871577419875E-5, 3.33373484773022E-5, 3.39130404042858E-6, 6.27885914770521E-5, 3.1309875198729E-5, 
                1.97209699963234E-5 )
    
        
#    deming.td9 <- pbreg(Y.td9~X.td9,data.frame(X.td9=X.td9,Y.td9=Y.td9), method=3)
    PBequi.small.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=F)
    PBequi.large.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=T)
        
    checkEquals( -2.44854969408334e-07, PBequi.small.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.small.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -2.44854969408334e-07, PBequi.large.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.large.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places

    # Theil-Sen
#    ts.td9 <- theilsen(Y.td9~X.td9,data.frame(X.td9=X.td9,Y.td9=Y.td9))
    ts.small.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=F,method.reg="TS")
    ts.large.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=T,method.reg="TS")
        
    checkEquals( 6.2991395442635e-07, ts.small.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.small.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 6.2991395442635e-07, ts.large.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.large.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    # TD13.xls data
    
    X.td13 <- c( 24794529.2480324, 54492735.891554, 71400324.331745, 8244180.28768875, 52509493.8140207, 49607163.0752203, 83827574.9515994, 47995388.3994996, 
                65394597.5719383, 60621910.5809957, 5483014.10574115, 84661259.5956848, 73871417.6242463, 47943755.9865663, 59617365.1220441, 59751718.1536454, 
                67526202.3095708, 43799674.6882773, 45633493.334196, 26732793.3066674, 53942813.5962028, 78076962.3282439, 96898323.9438709, 77073719.7541118, 
                29186602.9770011, 39530221.4111076, 78230070.4374343, 68911992.683726, 13516437.3215798, 28244384.8442956, 46735070.8937759, 79787020.4301856, 
                26196348.8304632, 67274300.1782288, 59058632.8684933, 34347561.2827838, 70726905.9209775, 13887668.6310818, 15243794.7219811, 51762357.5861341, 
                3777530.27538479, 44537496.2805357, 73284747.8866719, 24185258.6285496, 14458656.8718652, 31552195.7485487, 4041571.59760381, 72470788.8984188, 
                28884896.9586527, 18145471.5994012 )
        
    Y.td13 <- c( 21105237.1498578, 42316902.0455709, 65208707.2020304, 9247931.74778137, 66634360.9634913, 37957698.7793361, 95608490.3375025, 57881336.7248612, 
                59825287.298063, 62919716.7347258, 6453173.26688807, 68154558.3186772, 61432886.4897738, 36247951.2864404, 45366670.6876532, 45133759.6066822, 
                58880154.1326545, 53586686.3682414, 51087248.3563812, 20047078.7909169, 55296719.8327848, 73453100.6344544, 85666146.0292695, 69278207.448691, 
                27769147.8340553, 45079260.0802043, 74091530.765788, 85172813.2968523, 14816266.0519741, 28142372.2618154, 47079841.1549798, 92681460.6305504, 
                27362811.205277, 67593588.7363532, 59123504.6402203, 42074407.0566868, 58253866.2865107, 11517688.8782602, 13779895.5448975, 59500057.0904813, 
                3872350.51249896, 35939506.0007184, 100595978.033436, 19394350.6560089, 12787157.7419875, 33337348.4773022, 3391304.04042858, 62788591.4770521, 
                31309875.198729, 19720969.9963234 )
     
#     deming.td13 <- pbreg(Y.td13~X.td13,data.frame(X.td13=X.td13,Y.td13=Y.td13), method=3)
    PBequi.small.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=F)
    PBequi.large.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=T)
        
    checkEquals( -244854.969408333, PBequi.small.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.small.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -244854.969408333, PBequi.large.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.large.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
 
    #Theil-Sen
#    ts.td13 <- theilsen(Y.td13~X.td13,data.frame(X.td13=X.td13,Y.td13=Y.td13))
    ts.small.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=F,method.reg="TS")
    ts.large.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=T,method.reg="TS")
        
    checkEquals( 629913.954426358, ts.small.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.small.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 629913.954426358, ts.large.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.large.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
   
    
    
    # TD15.xls
    
    X.td15 <- c(18.654599109456, 17.0610823229594, 37.8076135618173, 51.1413961380522, 103.985667589592, 34.8605133072495, 35.1146708951621, 16.005996343635, 
                68.5890100912778, 62.1699013684068, 18.1641040703209, 94.2563842860796, 31.67278631354, 18.8059773112482, 67.9347403613216, 18.260239156489, 
                87.1550332248675, 98.5843421029197, 20.3747255514258, 4.20199417418543, 27.7662087802102, 16.6797331091928, 52.3335230459024, 1.25970439645986, 
                2.55285163008995, 71.1422300479687, 99.3428902581639, 40.6747911609727, 82.5924169862313, 64.7062991090824, 1.46968233573658, 83.5389530283147, 
                19.8368624090347, 44.339918470983, 19.2879672940529, 86.5484515948081, 79.8461255550168, 66.1164680506251, 32.8785828616434, 85.634465208222, 
                15.772867025856, 102.866851788539, 16.4442973704484, 23.0916021823352, 83.3011761592876, 9.47357732139885, 15.4574423263928, 67.0394769874472, 
                32.7239965276836, 71.7361088013058 )
        
    Y.td15 <- c(15.9116325401596, 17.9423397346645, 58.1755698848695, 49.0779434653576, 113.758268975683, 35.5525894416901, 39.4088027819502, 15.0167183931276, 
                65.8450151148211, 62.6667874119841, 15.8020862675435, 83.9418317622942, 33.5524291118479, 18.5775089262923, 73.872592767028, 21.600509295785, 
                93.4269960011703, 89.8278248717022, 20.130360709819, 4.94910345310486, 23.1566064984898, 17.336339530418, 56.2280624424048, 1.70969794990576, 
                3.40128933880261, 87.5197563754811, 67.3466761237754, 35.6882873995953, 90.4041198524321, 62.2603160684644, 1.52718831259026, 84.5412183739216, 
                19.0003036715583, 51.0624306952549, 12.0291575952956, 55.3970074896586, 96.7628055558332, 69.6940969840453, 29.1380990869214, 85.4597993067675, 
                18.1508502210932, 92.1385811560549, 17.012175857013, 23.0242126162252, 95.7510267835212, 8.70125548126123, 13.4955271740672, 45.8284262464427, 
                31.1175741976669, 69.5490311988868)
        
#     deming.td15 <- pbreg(Y.td15~X.td15,data.frame(X.td15=X.td15,Y.td15=Y.td15), method=3)
    PBequi.small.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=F)
    PBequi.large.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=T)
        
    checkEquals( -0.639752984809869, PBequi.small.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 1.02058852149448, PBequi.small.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.639752984809869, PBequi.large.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 1.02058852149448, PBequi.large.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#    ts.td15 <- theilsen(Y.td15~X.td15,data.frame(X.td15=X.td15,Y.td15=Y.td15))
    ts.small.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=F,method.reg="TS")
    ts.large.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=T,method.reg="TS")
        
    checkEquals( 0.215941922555261, ts.small.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.993801492789635, ts.small.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.215941922555261, ts.large.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.993801492789635, ts.large.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
   

    
    # TD20.xls
 
    X.td20 <- c(40.7416675085396, 83.2921878022294, 36.4584614937195, 72.0877727140495, 105.023225483845, 16.0054850259489, 77.939023513233, 67.9581547886773, 
                72.6256313765255, 76.1721718038233, 53.5660975152265, 44.5346723059445, 65.5050371740949, 34.7466804639362, 96.0978992891711, 22.2521041428996, 
                22.5836249521736, 62.0494626614161, 61.9702908007929, 28.2088726776012, 31.6304427667499, 44.0148198828374, 104.670709377516, 52.7016282328276, 
                1.60206831757677, 75.4250970942993, 9.34446754429201, 60.6762586794768, 73.933484555705, 92.9007223924939, 92.1358694560978, 17.1124285739673, 
                16.8171641673365, 54.6806652266699, 66.5315262775023, 92.626350509732, 89.9333484415825, 92.3201120961047, 76.0988174908423, 17.1096655807122, 
                67.8811974065475, 42.3956783976743, 40.0315902646466, 67.1597134632739, 50.835099276364, 35.8443954390043, 1.26426041082727, 41.0671736784787, 
                11.2645309896417, 99.9805522966942 )
     
     Y.td20 <- c(10.7521704713994, 13.4509410376188, 11.2856306611224, 10.7072031331434, 6.54511523826305, 6.19341973153328, 8.9913906536256, 12.3132146064145, 
                 13.0068162195109, 4.74699049037311, 9.23390774820934, 9.25776029197391, 9.76777184256182, 10.3206462496027, 7.63552298198159, 10.3097082932773, 
                 8.98004925615561, 13.0747427217869, 11.043558452402, 10.6615472307638, 9.06764709725244, 12.9103504115876, 10.9096799031437, 7.52693609703188, 
                 9.80669510453226, 10.0521231199264, 9.63095814224613, 8.5641931365438, 9.2882644554061, 10.6419750816148, 11.4130646454215, 11.4512769554377, 
                 13.2649021988222, 13.0654807233486, 8.23186894210709, 12.778014425933, 9.30152382689763, 9.08142311092212, 12.3867931320649, 11.1263164725831, 
                 7.55471439942213, 12.7413975083685, 8.94432971745523, 10.6831718241507, 9.95511449623938, 8.68553553767885, 8.5049120707315, 13.349972370573, 
                 10.3520245838025, 8.71363686499155)

    # pbreg from deming package only returns positive slope estimates, so change sign of Y!
#   deming.td20 <- pbreg(Y.td20~X.td20,data.frame(X.td20=X.td20,Y.td20=-Y.td20), method=3)

    PBequi.small.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=F)
    PBequi.large.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=T)
        
    checkEquals( 13.7983056265337, PBequi.small.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0654473128164153, PBequi.small.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 13.7983056265337, PBequi.large.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0654473128164153, PBequi.large.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#k    ts.td20 <- theilsen(Y.td20~X.td20,data.frame(X.td20=X.td20,Y.td20=Y.td20))
    ts.small.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=F,method.reg="TS")
    ts.large.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=T,method.reg="TS")
        
    checkEquals( 10.308472354129, ts.small.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0026117998040781, ts.small.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 10.308472354129, ts.large.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0026117998040781, ts.large.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    
    # TC41.xls
    X.td41 <- c(51.2645749891018, 25.057888623163, 87.7627598051843, 62.7082431751683, 80.1505333459012, 4.62656604009068, 33.8310555651344, 32.0639233133885, 
                55.2077551040875, 88.0356768073631, 26.4910934056359, 46.2676105006161, 67.5309426729014, 69.4845641448515, 99.6168230345125, 29.8740565635394, 
                76.4074905899151, 40.4037427291051, 102.582022432525, 77.2420615775276, 54.4091160177792, 31.6670196474807, 35.3688444029188, 48.4596331809692, 
                27.2766177138551, 16.8936814075242, 21.978946351695, 31.8737570958866, 65.2811331882764, 11.4235216747994, 26.5759954950747, 89.7483619905396, 
                70.8540190812346, 81.722195686426, 88.8717635815042, 69.8842397333777, 41.3339139273995, 46.7715712147959, 7.2758706475883, 107.885815963472, 
                72.476537675172, 5.47829692782881, 54.7122861461908, 100.774725860614, 58.9921395982599, 26.0127119387814, 52.4842022010358, 20.4930187951634, 
                19.2251154457486, 15.6062233212609)
        
    Y.td41 <- c(57.5184652931192, 20.8590706998877, 108.803321969406, 61.9186963049216, 66.3660888378279, 4.45211496523264, 32.5230480583973, 35.6080901524917, 
                53.8918711450082, 60.4588280634825, 26.36826855427, 41.3684710045443, 64.2210938093771, 56.4507675239752, 104.375984642575, 27.1349756417984, 
                94.57022726837, 29.4173054011996, 65.8469193441376, 72.8808353078608, 52.3451209927057, 29.9746067384427, 41.473958110454, 42.6744860392878, 
                36.4560529958421, 16.3614420035519, 25.0786009186081, 23.9554024584609, 71.1179574648198, 10.0486986665957, 19.0607414061787, 71.7842138052182, 
                126.559037624235, 94.0881361451604, 87.2145942516959, 89.5887499202462, 48.5784251203412, 35.6846235173564, 7.12184193955352, 100.323182519431, 
                69.23017359012, 3.82672858010298, 49.4382345176691, 77.5738341281733, 52.4666442933717, 14.6409084616032, 49.8976536701488, 21.1546025964904, 
                12.5229666712546, 13.2198613541858)
        
#    deming.td41 <- pbreg(Y.td41~X.td41,data.frame(X.td41=X.td41,Y.td41=Y.td41), method=3)
    PBequi.small.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=F)
    PBequi.large.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=T)
        
    checkEquals( -1.74861851923587, PBequi.small.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.996354386947289, PBequi.small.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -1.74861851923587, PBequi.large.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.996354386947289, PBequi.large.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#    ts.td41 <- theilsen(Y.td41~X.td41,data.frame(X.td41=X.td41,Y.td41=Y.td41))
    ts.small.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=F,method.reg="TS")
    ts.large.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=T,method.reg="TS")
        
    checkEquals( -0.33225717077892, ts.small.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.957048192302262, ts.small.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.33225717077892, ts.large.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.957048192302262, ts.large.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    # TD75.xls

    X.td75 <- c(90.9485400480344, 90.7431136750944, 26.3102304346173, 98.1923599874708, 35.0824973712106, 36.9969067737957, 14.8179253110874, 8.76477378922007, 
                32.2168398936375, 91.7976025914269, 75.2891547745564, 40.642084097967, 93.23660592041, 100.847125571886, 80.0367746433434, 9.73072817138115, 
                64.569921021703, 82.102749771495, 61.637913139448, 43.1726018718829, 30.9426405852636, 61.6489760843261, 72.5337200551608, 11.0971091511447, 
                91.1682684359647, 60.1522187636398, 53.6196723233099, 80.3991970785158, 46.7073650320478, 52.371013380707, 98.5451872718003, 31.1406239614125, 
                90.240170196133, 49.0589554374691, 86.9003789837115, 32.9579886227007, 24.0949524116379, 7.41009094921814, 81.1157940032336, 85.4694193753626, 
                42.1429310664425, 38.3807636449311, 19.8071119844752, 4.72048213046477, 86.48263405132, 88.9649734429405, 51.4221170355993, 88.8534761828361, 
                78.6355105622646, 11.799392585867)
    
    Y.td75 <- c(22.022678199682, 10.5849687064612, 70.545127040311, 1.91497891273671, 65.9280486240602, 70.4157396621642, 68.4545122022955, 91.022279815723, 
                52.5677072283795, 13.1589064936892, 29.1447890966354, 50.8408917564495, 19.7150011527249, 1.04103018782965, 17.6127620440527, 82.1291792966766, 
                42.250983349452, 5.60738627168227, 31.8795662581425, 54.3585165737907, 80.9769303190462, 40.7327484396029, 31.3717480824724, 94.1612819745663, 
                2.88484754250793, 46.6750774514867, 41.4465106338786, 20.031670042957, 64.765728125989, 51.7559200276114, 6.4110959903549, 44.664188884282, 
                6.12057655318134, 42.372324659659, 12.1204881694834, 65.8043448449057, 82.2468633319095, 104.200093777095, 22.7130359046475, 7.82813668061742, 
                53.9881500506756, 69.0313570401099, 67.425743171629, 92.1552906396942, 8.5364694548145, 17.1251097466213, 60.1446642607691, 18.3849543373421, 
                10.7104421280535, 83.3322475975441)

     # pbreg from deming package only returns positive slope estimates, so change sign of Y!      
#    deming.td75 <- pbreg(Y.td75~X.td75,data.frame(X.td75=X.td75,Y.td75=-Y.td75), method=3)
    PBequi.small.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=F)
    PBequi.large.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=T)
        
    checkEquals( 99.8952889829088, PBequi.small.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.995207226008001, PBequi.small.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 99.8952889829088, PBequi.large.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.995207226008001, PBequi.large.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places

    #Theil-Sen
#    ts.td75 <- theilsen(Y.td75~X.td75,data.frame(X.td75=X.td75,Y.td75=Y.td75))
    ts.small.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=F,method.reg="TS")
    ts.large.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=T,method.reg="TS")
        
    checkEquals( 96.629630145324, ts.small.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.947856465074523, ts.small.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 96.629630145324, ts.large.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.947856465074523, ts.large.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    
}

test.mc.PBequi.resampling1 <- function()
{
	
	# Check problematic uncorrelated data set for borderline case where shift index K is half the sample size
	# this would result in index error with original PaBa equations. Should be properly handled by our PaBa implementation.
	
	td.idxprob <- matrix(ncol=2,byrow=T,data=c(1.25,1.81,1.11,1.15,1.67,0.78,1.22,1.53,1.33,1.22,0.99,1.06,1.40,1.30,1.34,0.82,0.53,2.30,1.52,2.11))
	
	set.seed(42)
	f <- function(x,y) {
		ind <- sample(1:length(x), length(x), replace=TRUE)
		test1 <- try( mcr:::mc.PBequi(x[ind],y[ind], calcCI=FALSE) )		# test without CI
		checkTrue( class(test1) != "try-error")
		test2 <- try( mcr:::mc.PBequi(x[ind],y[ind], calcCI=TRUE) )		# test with CI
		checkTrue( class(test2) != "try-error")
	}
	
	for(i in 1:2000) {f(td.idxprob[,1],td.idxprob[,2])} # call mc.PBequi with different data configurations, first occurence of K=N/2 configuration in iteration 734


	# Check another problematic data constellation with seed forcing an error in the implementation untill version <=V1.2.
	# The problem occurs due to differences in the correlation of global data and resampled data when the angle matrix
	# is resampled from global angle matrix but correlation of method changes in sign.
	
	# run via interface
	Xval <- c(4, 4.7, 6.1, 5.8, 7.2, 1.9, 4.2, 5.2, 4.8, 5.1, 5.5, 4.6, 5.4, 4.2)
	Yval <- c(4.29, 5.08, 1.82, 2.12, 2.21, 2.01, 4.51, 1.44, 5.05, 5.53, 5.77, 1.64, 5.66, 4.52)
	

	fit1 <- try( mcreg(Xval, Yval, method.reg="PBequi", method.ci="bootstrap", method.bootstrap.ci="quantile", rng.seed=1458260) )	
	checkTrue(class(fit1) != "try-error")				# no error occured	
}



cat("\n\nmcPBequi.R method comparison test cases\n\n")

## check values for equiariant Passing-Bablok and Theil-Sen regression against regression coefficients computed with the deming package (v. 1.4) by Terry Therneau. 
## Confidence intervals are calculated by another method and are therefore not directly comparable. 

#library(deming)

test.mc.PBequi.call2 <- function() 
{
    data(creatinine)
    crea <- creatinine[complete.cases(creatinine),]

# Compare to deming package
#    library(deming)

#    Compare equivariant PaBa regression

#    pbdeming <- pbreg(serum.crea~plasma.crea,crea, method=3)
#    dput(pbdeming[[1]])
#Output: c(`(Intercept)` = 0.102307692307692, plasma.crea = 0.923076923076924 

    pbeqsmall <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="PBequi",calcCI=TRUE,methodlarge=F)
    pbeqlarge <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="PBequi",calcCI=TRUE,methodlarge=T)
    
    
    
    
    # values are copied and pasted from deming package
    
                   
    checkEquals( 0.102307692307692, pbeqsmall$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.923076923076924, pbeqsmall$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 14 decimal places
     checkEquals( 0.102307692307692, pbeqlarge$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 14 decimal places
    checkEquals( 0.923076923076924, pbeqlarge$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 14 decimal places

# Compare Theil-Sen regression

#     tsdeming <- theilsen(serum.crea~plasma.crea,crea )  
#     dput(tsdeming[[1]])   
     tssmall <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="TS",calcCI=TRUE,methodlarge=F)
    tslarge <- mcr:::mc.PBequi(crea$plasma.crea,crea$serum.crea,slope.measure="radian",method.reg="TS",calcCI=TRUE,methodlarge=T)

    checkEquals( 0.17796875, tssmall$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.84375, tssmall$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
     checkEquals( 0.17796875, tslarge$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.84375, tslarge$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
  
    X.td9 <- c( 2.47945292480324E-5, 5.4492735891554E-5, 7.1400324331745E-5, 8.24418028768875E-6, 5.25094938140207E-5, 4.96071630752203E-5, 8.38275749515994E-5, 
                4.79953883994996E-5, 6.53945975719383E-5, 6.06219105809957E-5, 5.48301410574115E-6, 8.46612595956848E-5, 7.38714176242463E-5, 4.79437559865663E-5, 
                5.96173651220441E-5, 5.97517181536454E-5, 6.75262023095708E-5, 4.37996746882773E-5, 4.5633493334196E-5, 2.67327933066674E-5, 5.39428135962028E-5, 
                7.80769623282439E-5, 9.68983239438709E-5, 7.70737197541118E-5, 2.91866029770011E-5, 3.95302214111076E-5, 7.82300704374343E-5, 6.8911992683726E-5, 
                1.35164373215798E-5, 2.82443848442956E-5, 4.67350708937759E-5, 7.97870204301856E-5, 2.61963488304632E-5, 6.72743001782288E-5, 5.90586328684933E-5, 
                3.43475612827838E-5, 7.07269059209775E-5, 1.38876686310818E-5, 1.52437947219811E-5, 5.17623575861341E-5, 3.77753027538479E-6, 4.45374962805357E-5, 
                7.32847478866719E-5, 2.41852586285496E-5, 1.44586568718652E-5, 3.15521957485487E-5, 4.04157159760381E-6, 7.24707888984188E-5, 2.88848969586527E-5, 
                1.81454715994012E-5 )
    
    Y.td9 <- c( 2.11052371498578E-5, 4.23169020455709E-5, 6.52087072020304E-5, 9.24793174778137E-6, 6.66343609634913E-5, 3.79576987793361E-5, 9.56084903375025E-5, 
                5.78813367248612E-5, 5.9825287298063E-5, 6.29197167347258E-5, 6.45317326688806E-6, 6.81545583186772E-5, 6.14328864897738E-5, 3.62479512864404E-5, 
                4.53666706876532E-5, 4.51337596066822E-5, 5.88801541326545E-5, 5.35866863682414E-5, 5.10872483563812E-5, 2.00470787909169E-5, 5.52967198327848E-5, 
                7.34531006344544E-5, 8.56661460292695E-5, 6.9278207448691E-5, 2.77691478340553E-5, 4.50792600802043E-5, 7.4091530765788E-5, 8.51728132968523E-5, 
                1.48162660519741E-5, 2.81423722618154E-5, 4.70798411549798E-5, 9.26814606305503E-5, 2.7362811205277E-5, 6.75935887363532E-5, 5.91235046402203E-5, 
                4.20744070566868E-5, 5.82538662865107E-5, 1.15176888782602E-5, 1.37798955448975E-5, 5.95000570904813E-5, 3.87235051249896E-6, 3.59395060007184E-5, 
                0.000100595978033436, 1.93943506560089E-5, 1.27871577419875E-5, 3.33373484773022E-5, 3.39130404042858E-6, 6.27885914770521E-5, 3.1309875198729E-5, 
                1.97209699963234E-5 )
    
        
#    deming.td9 <- pbreg(Y.td9~X.td9,data.frame(X.td9=X.td9,Y.td9=Y.td9), method=3)
    PBequi.small.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=F)
    PBequi.large.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=T)
        
    checkEquals( -2.44854969408334e-07, PBequi.small.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.small.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -2.44854969408334e-07, PBequi.large.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.large.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places

    # Theil-Sen
#    ts.td9 <- theilsen(Y.td9~X.td9,data.frame(X.td9=X.td9,Y.td9=Y.td9))
    ts.small.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=F,method.reg="TS")
    ts.large.td9 <- mcr:::mc.PBequi(X.td9, Y.td9,methodlarge=T,method.reg="TS")
        
    checkEquals( 6.2991395442635e-07, ts.small.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.small.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 6.2991395442635e-07, ts.large.td9$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.large.td9$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    # TD13.xls data
    
    X.td13 <- c( 24794529.2480324, 54492735.891554, 71400324.331745, 8244180.28768875, 52509493.8140207, 49607163.0752203, 83827574.9515994, 47995388.3994996, 
                65394597.5719383, 60621910.5809957, 5483014.10574115, 84661259.5956848, 73871417.6242463, 47943755.9865663, 59617365.1220441, 59751718.1536454, 
                67526202.3095708, 43799674.6882773, 45633493.334196, 26732793.3066674, 53942813.5962028, 78076962.3282439, 96898323.9438709, 77073719.7541118, 
                29186602.9770011, 39530221.4111076, 78230070.4374343, 68911992.683726, 13516437.3215798, 28244384.8442956, 46735070.8937759, 79787020.4301856, 
                26196348.8304632, 67274300.1782288, 59058632.8684933, 34347561.2827838, 70726905.9209775, 13887668.6310818, 15243794.7219811, 51762357.5861341, 
                3777530.27538479, 44537496.2805357, 73284747.8866719, 24185258.6285496, 14458656.8718652, 31552195.7485487, 4041571.59760381, 72470788.8984188, 
                28884896.9586527, 18145471.5994012 )
        
    Y.td13 <- c( 21105237.1498578, 42316902.0455709, 65208707.2020304, 9247931.74778137, 66634360.9634913, 37957698.7793361, 95608490.3375025, 57881336.7248612, 
                59825287.298063, 62919716.7347258, 6453173.26688807, 68154558.3186772, 61432886.4897738, 36247951.2864404, 45366670.6876532, 45133759.6066822, 
                58880154.1326545, 53586686.3682414, 51087248.3563812, 20047078.7909169, 55296719.8327848, 73453100.6344544, 85666146.0292695, 69278207.448691, 
                27769147.8340553, 45079260.0802043, 74091530.765788, 85172813.2968523, 14816266.0519741, 28142372.2618154, 47079841.1549798, 92681460.6305504, 
                27362811.205277, 67593588.7363532, 59123504.6402203, 42074407.0566868, 58253866.2865107, 11517688.8782602, 13779895.5448975, 59500057.0904813, 
                3872350.51249896, 35939506.0007184, 100595978.033436, 19394350.6560089, 12787157.7419875, 33337348.4773022, 3391304.04042858, 62788591.4770521, 
                31309875.198729, 19720969.9963234 )
     
#     deming.td13 <- pbreg(Y.td13~X.td13,data.frame(X.td13=X.td13,Y.td13=Y.td13), method=3)
    PBequi.small.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=F)
    PBequi.large.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=T)
        
    checkEquals( -244854.969408333, PBequi.small.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.small.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -244854.969408333, PBequi.large.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.991592719686709, PBequi.large.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
 
    #Theil-Sen
#    ts.td13 <- theilsen(Y.td13~X.td13,data.frame(X.td13=X.td13,Y.td13=Y.td13))
    ts.small.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=F,method.reg="TS")
    ts.large.td13 <- mcr:::mc.PBequi(X.td13, Y.td13,methodlarge=T,method.reg="TS")
        
    checkEquals( 629913.954426358, ts.small.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.small.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 629913.954426358, ts.large.td13$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.936547814163375, ts.large.td13$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
   
    
    
    # TD15.xls
    
    X.td15 <- c(18.654599109456, 17.0610823229594, 37.8076135618173, 51.1413961380522, 103.985667589592, 34.8605133072495, 35.1146708951621, 16.005996343635, 
                68.5890100912778, 62.1699013684068, 18.1641040703209, 94.2563842860796, 31.67278631354, 18.8059773112482, 67.9347403613216, 18.260239156489, 
                87.1550332248675, 98.5843421029197, 20.3747255514258, 4.20199417418543, 27.7662087802102, 16.6797331091928, 52.3335230459024, 1.25970439645986, 
                2.55285163008995, 71.1422300479687, 99.3428902581639, 40.6747911609727, 82.5924169862313, 64.7062991090824, 1.46968233573658, 83.5389530283147, 
                19.8368624090347, 44.339918470983, 19.2879672940529, 86.5484515948081, 79.8461255550168, 66.1164680506251, 32.8785828616434, 85.634465208222, 
                15.772867025856, 102.866851788539, 16.4442973704484, 23.0916021823352, 83.3011761592876, 9.47357732139885, 15.4574423263928, 67.0394769874472, 
                32.7239965276836, 71.7361088013058 )
        
    Y.td15 <- c(15.9116325401596, 17.9423397346645, 58.1755698848695, 49.0779434653576, 113.758268975683, 35.5525894416901, 39.4088027819502, 15.0167183931276, 
                65.8450151148211, 62.6667874119841, 15.8020862675435, 83.9418317622942, 33.5524291118479, 18.5775089262923, 73.872592767028, 21.600509295785, 
                93.4269960011703, 89.8278248717022, 20.130360709819, 4.94910345310486, 23.1566064984898, 17.336339530418, 56.2280624424048, 1.70969794990576, 
                3.40128933880261, 87.5197563754811, 67.3466761237754, 35.6882873995953, 90.4041198524321, 62.2603160684644, 1.52718831259026, 84.5412183739216, 
                19.0003036715583, 51.0624306952549, 12.0291575952956, 55.3970074896586, 96.7628055558332, 69.6940969840453, 29.1380990869214, 85.4597993067675, 
                18.1508502210932, 92.1385811560549, 17.012175857013, 23.0242126162252, 95.7510267835212, 8.70125548126123, 13.4955271740672, 45.8284262464427, 
                31.1175741976669, 69.5490311988868)
        
#     deming.td15 <- pbreg(Y.td15~X.td15,data.frame(X.td15=X.td15,Y.td15=Y.td15), method=3)
    PBequi.small.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=F)
    PBequi.large.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=T)
        
    checkEquals( -0.639752984809869, PBequi.small.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 1.02058852149448, PBequi.small.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.639752984809869, PBequi.large.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 1.02058852149448, PBequi.large.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#    ts.td15 <- theilsen(Y.td15~X.td15,data.frame(X.td15=X.td15,Y.td15=Y.td15))
    ts.small.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=F,method.reg="TS")
    ts.large.td15 <- mcr:::mc.PBequi(X.td15, Y.td15,methodlarge=T,method.reg="TS")
        
    checkEquals( 0.215941922555261, ts.small.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.993801492789635, ts.small.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.215941922555261, ts.large.td15$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.993801492789635, ts.large.td15$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
   

    
    # TD20.xls
 
    X.td20 <- c(40.7416675085396, 83.2921878022294, 36.4584614937195, 72.0877727140495, 105.023225483845, 16.0054850259489, 77.939023513233, 67.9581547886773, 
                72.6256313765255, 76.1721718038233, 53.5660975152265, 44.5346723059445, 65.5050371740949, 34.7466804639362, 96.0978992891711, 22.2521041428996, 
                22.5836249521736, 62.0494626614161, 61.9702908007929, 28.2088726776012, 31.6304427667499, 44.0148198828374, 104.670709377516, 52.7016282328276, 
                1.60206831757677, 75.4250970942993, 9.34446754429201, 60.6762586794768, 73.933484555705, 92.9007223924939, 92.1358694560978, 17.1124285739673, 
                16.8171641673365, 54.6806652266699, 66.5315262775023, 92.626350509732, 89.9333484415825, 92.3201120961047, 76.0988174908423, 17.1096655807122, 
                67.8811974065475, 42.3956783976743, 40.0315902646466, 67.1597134632739, 50.835099276364, 35.8443954390043, 1.26426041082727, 41.0671736784787, 
                11.2645309896417, 99.9805522966942 )
     
     Y.td20 <- c(10.7521704713994, 13.4509410376188, 11.2856306611224, 10.7072031331434, 6.54511523826305, 6.19341973153328, 8.9913906536256, 12.3132146064145, 
                 13.0068162195109, 4.74699049037311, 9.23390774820934, 9.25776029197391, 9.76777184256182, 10.3206462496027, 7.63552298198159, 10.3097082932773, 
                 8.98004925615561, 13.0747427217869, 11.043558452402, 10.6615472307638, 9.06764709725244, 12.9103504115876, 10.9096799031437, 7.52693609703188, 
                 9.80669510453226, 10.0521231199264, 9.63095814224613, 8.5641931365438, 9.2882644554061, 10.6419750816148, 11.4130646454215, 11.4512769554377, 
                 13.2649021988222, 13.0654807233486, 8.23186894210709, 12.778014425933, 9.30152382689763, 9.08142311092212, 12.3867931320649, 11.1263164725831, 
                 7.55471439942213, 12.7413975083685, 8.94432971745523, 10.6831718241507, 9.95511449623938, 8.68553553767885, 8.5049120707315, 13.349972370573, 
                 10.3520245838025, 8.71363686499155)

    # pbreg from deming package only returns positive slope estimates, so change sign of Y!
#   deming.td20 <- pbreg(Y.td20~X.td20,data.frame(X.td20=X.td20,Y.td20=-Y.td20), method=3)

    PBequi.small.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=F)
    PBequi.large.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=T)
        
    checkEquals( 13.7983056265337, PBequi.small.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0654473128164153, PBequi.small.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 13.7983056265337, PBequi.large.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0654473128164153, PBequi.large.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#k    ts.td20 <- theilsen(Y.td20~X.td20,data.frame(X.td20=X.td20,Y.td20=Y.td20))
    ts.small.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=F,method.reg="TS")
    ts.large.td20 <- mcr:::mc.PBequi(X.td20, Y.td20,methodlarge=T,method.reg="TS")
        
    checkEquals( 10.308472354129, ts.small.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0026117998040781, ts.small.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 10.308472354129, ts.large.td20$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.0026117998040781, ts.large.td20$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    
    # TC41.xls
    X.td41 <- c(51.2645749891018, 25.057888623163, 87.7627598051843, 62.7082431751683, 80.1505333459012, 4.62656604009068, 33.8310555651344, 32.0639233133885, 
                55.2077551040875, 88.0356768073631, 26.4910934056359, 46.2676105006161, 67.5309426729014, 69.4845641448515, 99.6168230345125, 29.8740565635394, 
                76.4074905899151, 40.4037427291051, 102.582022432525, 77.2420615775276, 54.4091160177792, 31.6670196474807, 35.3688444029188, 48.4596331809692, 
                27.2766177138551, 16.8936814075242, 21.978946351695, 31.8737570958866, 65.2811331882764, 11.4235216747994, 26.5759954950747, 89.7483619905396, 
                70.8540190812346, 81.722195686426, 88.8717635815042, 69.8842397333777, 41.3339139273995, 46.7715712147959, 7.2758706475883, 107.885815963472, 
                72.476537675172, 5.47829692782881, 54.7122861461908, 100.774725860614, 58.9921395982599, 26.0127119387814, 52.4842022010358, 20.4930187951634, 
                19.2251154457486, 15.6062233212609)
        
    Y.td41 <- c(57.5184652931192, 20.8590706998877, 108.803321969406, 61.9186963049216, 66.3660888378279, 4.45211496523264, 32.5230480583973, 35.6080901524917, 
                53.8918711450082, 60.4588280634825, 26.36826855427, 41.3684710045443, 64.2210938093771, 56.4507675239752, 104.375984642575, 27.1349756417984, 
                94.57022726837, 29.4173054011996, 65.8469193441376, 72.8808353078608, 52.3451209927057, 29.9746067384427, 41.473958110454, 42.6744860392878, 
                36.4560529958421, 16.3614420035519, 25.0786009186081, 23.9554024584609, 71.1179574648198, 10.0486986665957, 19.0607414061787, 71.7842138052182, 
                126.559037624235, 94.0881361451604, 87.2145942516959, 89.5887499202462, 48.5784251203412, 35.6846235173564, 7.12184193955352, 100.323182519431, 
                69.23017359012, 3.82672858010298, 49.4382345176691, 77.5738341281733, 52.4666442933717, 14.6409084616032, 49.8976536701488, 21.1546025964904, 
                12.5229666712546, 13.2198613541858)
        
#    deming.td41 <- pbreg(Y.td41~X.td41,data.frame(X.td41=X.td41,Y.td41=Y.td41), method=3)
    PBequi.small.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=F)
    PBequi.large.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=T)
        
    checkEquals( -1.74861851923587, PBequi.small.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.996354386947289, PBequi.small.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -1.74861851923587, PBequi.large.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.996354386947289, PBequi.large.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    #Theil-Sen
#    ts.td41 <- theilsen(Y.td41~X.td41,data.frame(X.td41=X.td41,Y.td41=Y.td41))
    ts.small.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=F,method.reg="TS")
    ts.large.td41 <- mcr:::mc.PBequi(X.td41, Y.td41,methodlarge=T,method.reg="TS")
        
    checkEquals( -0.33225717077892, ts.small.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.957048192302262, ts.small.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.33225717077892, ts.large.td41$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 0.957048192302262, ts.large.td41$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    # TD75.xls

    X.td75 <- c(90.9485400480344, 90.7431136750944, 26.3102304346173, 98.1923599874708, 35.0824973712106, 36.9969067737957, 14.8179253110874, 8.76477378922007, 
                32.2168398936375, 91.7976025914269, 75.2891547745564, 40.642084097967, 93.23660592041, 100.847125571886, 80.0367746433434, 9.73072817138115, 
                64.569921021703, 82.102749771495, 61.637913139448, 43.1726018718829, 30.9426405852636, 61.6489760843261, 72.5337200551608, 11.0971091511447, 
                91.1682684359647, 60.1522187636398, 53.6196723233099, 80.3991970785158, 46.7073650320478, 52.371013380707, 98.5451872718003, 31.1406239614125, 
                90.240170196133, 49.0589554374691, 86.9003789837115, 32.9579886227007, 24.0949524116379, 7.41009094921814, 81.1157940032336, 85.4694193753626, 
                42.1429310664425, 38.3807636449311, 19.8071119844752, 4.72048213046477, 86.48263405132, 88.9649734429405, 51.4221170355993, 88.8534761828361, 
                78.6355105622646, 11.799392585867)
    
    Y.td75 <- c(22.022678199682, 10.5849687064612, 70.545127040311, 1.91497891273671, 65.9280486240602, 70.4157396621642, 68.4545122022955, 91.022279815723, 
                52.5677072283795, 13.1589064936892, 29.1447890966354, 50.8408917564495, 19.7150011527249, 1.04103018782965, 17.6127620440527, 82.1291792966766, 
                42.250983349452, 5.60738627168227, 31.8795662581425, 54.3585165737907, 80.9769303190462, 40.7327484396029, 31.3717480824724, 94.1612819745663, 
                2.88484754250793, 46.6750774514867, 41.4465106338786, 20.031670042957, 64.765728125989, 51.7559200276114, 6.4110959903549, 44.664188884282, 
                6.12057655318134, 42.372324659659, 12.1204881694834, 65.8043448449057, 82.2468633319095, 104.200093777095, 22.7130359046475, 7.82813668061742, 
                53.9881500506756, 69.0313570401099, 67.425743171629, 92.1552906396942, 8.5364694548145, 17.1251097466213, 60.1446642607691, 18.3849543373421, 
                10.7104421280535, 83.3322475975441)

     # pbreg from deming package only returns positive slope estimates, so change sign of Y!      
#    deming.td75 <- pbreg(Y.td75~X.td75,data.frame(X.td75=X.td75,Y.td75=-Y.td75), method=3)
    PBequi.small.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=F)
    PBequi.large.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=T)
        
    checkEquals( 99.8952889829088, PBequi.small.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.995207226008001, PBequi.small.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 99.8952889829088, PBequi.large.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.995207226008001, PBequi.large.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places

    #Theil-Sen
#    ts.td75 <- theilsen(Y.td75~X.td75,data.frame(X.td75=X.td75,Y.td75=Y.td75))
    ts.small.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=F,method.reg="TS")
    ts.large.td75 <- mcr:::mc.PBequi(X.td75, Y.td75,methodlarge=T,method.reg="TS")
        
    checkEquals( 96.629630145324, ts.small.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.947856465074523, ts.small.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( 96.629630145324, ts.large.td75$b0, tolerance=1e-12) 
    # Note: There is equivalence of results checked up to 12 decimal places
    checkEquals( -0.947856465074523, ts.large.td75$b1, tolerance=1e-12)
    # Note: There is equivalence of results checked up to 12 decimal places
    
    
    
}

test.mc.PBequi.resampling2 <- function()
{
	
	# Check problematic uncorrelated data set for borderline case where shift index K is half the sample size
	# this would result in index error with original PaBa equations. Should be properly handled by our PaBa implementation.
	
	td.idxprob <- matrix(ncol=2,byrow=T,data=c(1.25,1.81,1.11,1.15,1.67,0.78,1.22,1.53,1.33,1.22,0.99,1.06,1.40,1.30,1.34,0.82,0.53,2.30,1.52,2.11))
	
	set.seed(42)
	f <- function(x,y) {
		ind <- sample(1:length(x), length(x), replace=TRUE)
		test1 <- try( mcr:::mc.PBequi(x[ind],y[ind], calcCI=FALSE) )		# test without CI
		checkTrue( class(test1) != "try-error")
		test2 <- try( mcr:::mc.PBequi(x[ind],y[ind], calcCI=TRUE) )		# test with CI
		checkTrue( class(test2) != "try-error")
	}
	
	for(i in 1:2000) {f(td.idxprob[,1],td.idxprob[,2])} # call mc.PBequi with different data configurations, first occurence of K=N/2 configuration in iteration 734


	# Check another problematic data constellation with seed forcing an error in the implementation untill version <=V1.2.
	# The problem occurs due to differences in the correlation of global data and resampled data when the angle matrix
	# is resampled from global angle matrix but correlation of method changes in sign.
	
	# run via interface
	Xval <- c(4, 4.7, 6.1, 5.8, 7.2, 1.9, 4.2, 5.2, 4.8, 5.1, 5.5, 4.6, 5.4, 4.2)
	Yval <- c(4.29, 5.08, 1.82, 2.12, 2.21, 2.01, 4.51, 1.44, 5.05, 5.53, 5.77, 1.64, 5.66, 4.52)
	

	fit1 <- try( mcreg(Xval, Yval, method.reg="PBequi", method.ci="bootstrap", method.bootstrap.ci="quantile", rng.seed=1458260) )	
	checkTrue(class(fit1) != "try-error")				# no error occured	
}



test.mc.PBequi.pabasim <- function()
{

pabasim <- function(nsim=1, n=100, inter0 = 1, m0=2, sd0=0.1, pow=1, rmin=0, rmax=1, alpha=0.05,ylim=c(-0.05,0.05), slope.measure=c("radian","tangent"), method.reg=c("PBequi","TS"),methodlarge=F){
    # this function compares the means of the parameters obtained to the ones calculated from repeated simulations
    z <- qt(1-alpha/2,n-2)
	slope.measure <- match.arg(slope.measure)
	method.reg <- match.arg(method.reg)
    medint <- 1:nsim
    medintm <- 1:nsim
    tau <- 1:nsim
    zeta <- 1:nsim
    x0 <-1:nsim
    ncls <- 9
    results <- matrix(rep(0,nsim*ncls),ncol=ncls,nrow=nsim)

    for (i in 1:nsim){
        set.seed(i)
        #isim <<-i
        #print(paste0("simulation run: ",i))
        r <- (runif(n,rmin,rmax))
        if (method.reg=="PBequi"){
            X <- r+r^pow*rnorm(n,0,sd0)
        }else{
            X <- r
        }
        Y <- inter0+m0*(r+r^pow*rnorm(n,0,sd0))
        res <- mcr:::mc.PBequi(X,Y,method.reg=method.reg,slope.measure=slope.measure,extended.output=T,methodlarge=methodlarge,calcCI=T)
        results[i,] <- c(res[[1]],res[[2]],res[[3]],res[[4]],res[[5]],res[[6]],res[[7]],res[[8]],res[[9]])
        medint[i] <- median(Y-m0*X)
        if(method.reg=="PBequi"){
            tau[i] <- cor(Y-m0*X,Y+m0*X,method="kendall")
        }else{
            tau[i] <- cor(Y-m0*X,X,method="kendall")
        }
        zeta[i] <- sum(sign(Y-m0*X-inter0))
    }
#    # test intercept
#    print(paste0("True intercept: ",inter0, "; mean estimated intercept: ", mean(results[,1])))
#    #test slope
#    print(paste0("True slope: ",m0, "; mean estimated slope: ", mean(results[,2])))
#    #test sd of slope 
#    print(paste0("Empirical SD of slopes: ",sd(results[,2]), "; mean estimated sd: ", mean(results[,4])))
#    #method to determine the quantiles differs slightly between quadratic and quasilinear method!
#    # test sd of intercept for fixed slope m0 
#    print(paste0("Empirical SD of intercepts: ",sd(results[,1]), "; mean estimated sd: ", mean(results[,3])))
#    #method to determine the quantiles differs slightly between quadratic and quasilinear method. This propagates also into the estimation of the intercept sd. 
#    # test sd of zeta
#    print(paste0("Empirical SD of zeta: ",sd(zeta), "; mean estimated sd: ", sqrt(n)))
#    # test sd of tau
#    print(paste0("Empirical SD of tau: ",sd(tau), "; mean estimated sd: ", mean(sqrt(results[,7]))))
#    # test cor 
#    print(paste0("Empirical cor of tau and zeta: ",cor(tau,zeta,method="pearson"), "; mean estimated cor: ", mean(results[,8])))
#    #test sd of intercept for fixed m 
#    print(paste0("Empirical SD of intercept for fixed slope: ",sd(medint), "; mean estimated sd: ", mean(results[,6])))
    

    if(method.reg=="PBequi"){
        checkEquals( 1.00141900097656, mean(results[,1]), tolerance=1e-12) 
        checkEquals( 4.48916413425509, mean(results[,2]), tolerance=1e-12) 
        checkEquals( 0.117864362025096, sd(results[,2]), tolerance=1e-12) 
        if(methodlarge==F){
            checkEquals( 0.120603800379936, mean(results[,4]), tolerance=1e-12) 
        }else{
            checkEquals( 0.120613070650643, mean(results[,4]), tolerance=1e-12) 
        }
        checkEquals( 0.0174845691829908, sd(results[,1]), tolerance=1e-12) 
        if(methodlarge==F){
            checkEquals( 0.0257175181285312, mean(results[,3]), tolerance=1e-12) 
        }else{
            checkEquals( 0.0257165435814335, mean(results[,3]), tolerance=1e-12) 
        }
        checkEquals( 8.85061203156784, sd(zeta), tolerance=1e-12) 
        checkEquals( 0.0697719694242071, sd(tau), tolerance=1e-12) 
        checkEquals( 0.0712810095037541, mean(sqrt(results[,7])), tolerance=1e-12) 
        checkEquals( 0.687031297354726, cor(tau,zeta,method="pearson"), tolerance=1e-12) 
        checkEquals( 0.598614714856513, mean(results[,8]), tolerance=1e-12) 
        checkEquals( 0.018227606026415, sd(medint), tolerance=1e-12) 
        checkEquals( 0.0190851670031283, mean(results[,6]), tolerance=1e-12) 
    }else{
        checkEquals( 0.998277096719531, mean(results[,1]), tolerance=1e-12) 
        checkEquals( 4.51571430483745, mean(results[,2]), tolerance=1e-12) 
        checkEquals( 0.09160094461838, sd(results[,2]), tolerance=1e-12) 
        if(methodlarge==F){
            checkEquals( 0.0877858229514845, mean(results[,4]), tolerance=1e-12) 
        }else{
            checkEquals( 0.0877925005149715, mean(results[,4]), tolerance=1e-12) 
        }
        checkEquals( 0.0168999239835992, sd(results[,1]), tolerance=1e-12) 
        if(methodlarge==F){
            checkEquals( 0.0196454687476921, mean(results[,3]), tolerance=1e-12) 
        }else{
            checkEquals( 0.0196465215475306, mean(results[,3]), tolerance=1e-12) 
        }
        checkEquals( 10.430627008861, sd(zeta), tolerance=1e-12) 
        checkEquals( 0.0764477234760266, sd(tau), tolerance=1e-12) 
        checkEquals( 0.0717175331962926, mean(sqrt(results[,7])), tolerance=1e-12) 
        checkEquals( 0.690133407745362, cor(tau,zeta,method="pearson"), tolerance=1e-12) 
        checkEquals( 0.592228200366131, mean(results[,8]), tolerance=1e-12) 
        checkEquals( 0.0137260093560068, sd(medint), tolerance=1e-12) 
        checkEquals( 0.0137009842376418, mean(results[,6]), tolerance=1e-12) 
    }
}

pabasim(nsim=100,n=100,rmin=0,pow=1,m0=4.5,sd0=0.1, method.reg="PBequi",ylim=c(-0.1,0.1),slope.measure="tangent",methodlarge=F)
pabasim(nsim=100,n=100,rmin=0,pow=1,m0=4.5,sd0=0.1, method.reg="PBequi",ylim=c(-0.1,0.1),slope.measure="tangent",methodlarge=T)
pabasim(nsim=100,n=100,rmin=0,pow=1,m0=4.5,sd0=0.1, method.reg="TS",ylim=c(-0.1,0.1),slope.measure="tangent",methodlarge=F)
pabasim(nsim=100,n=100,rmin=0,pow=1,m0=4.5,sd0=0.1, method.reg="TS",ylim=c(-0.1,0.1),slope.measure="tangent",methodlarge=T)
}