R/epi.ltd.R

In epiR: Tools for the Analysis of Epidemiological Data

epi.ltd <- function(dat, std = "305"){

if (std == "305"){

std <- 305

a0 <- rep(1, times = 305)

a1 <- c(-1,-0.99342,-0.98684,-0.98026,-0.97368,-0.96711,-0.96053,-0.95395,-0.94737,
-0.94079,-0.93421,-0.92763,-0.92105,-0.91447,-0.90789,-0.90132,-0.89474,
-0.88816,-0.88158,-0.875,-0.86842,-0.86184,-0.85526,-0.84868,-0.84211,
-0.83553,-0.82895,-0.82237,-0.81579,-0.80921,-0.80263,-0.79605,-0.78947,
-0.78289,-0.77632,-0.76974,-0.76316,-0.75658,-0.75,-0.74342,-0.73684,
-0.73026,-0.72368,-0.71711,-0.71053,-0.70395,-0.69737,-0.69079,-0.68421,
-0.67763,-0.67105,-0.66447,-0.65789,-0.65132,-0.64474,-0.63816,-0.63158,
-0.625,-0.61842,-0.61184,-0.60526,-0.59868,-0.59211,-0.58553,-0.57895,
-0.57237,-0.56579,-0.55921,-0.55263,-0.54605,-0.53947,-0.53289,-0.52632,
-0.51974,-0.51316,-0.50658,-0.5,-0.49342,-0.48684,-0.48026,-0.47368,
-0.46711,-0.46053,-0.45395,-0.44737,-0.44079,-0.43421,-0.42763,-0.42105,
-0.41447,-0.40789,-0.40132,-0.39474,-0.38816,-0.38158,-0.375,-0.36842,
-0.36184,-0.35526,-0.34868,-0.34211,-0.33553,-0.32895,-0.32237,-0.31579,
-0.30921,-0.30263,-0.29605,-0.28947,-0.28289,-0.27632,-0.26974,-0.26316,
-0.25658,-0.25,-0.24342,-0.23684,-0.23026,-0.22368,-0.21711,-0.21053,
-0.20395,-0.19737,-0.19079,-0.18421,-0.17763,-0.17105,-0.16447,-0.15789,
-0.15132,-0.14474,-0.13816,-0.13158,-0.125,-0.11842,-0.11184,-0.10526,
-0.09868,-0.09211,-0.08553,-0.07895,-0.07237,-0.06579,-0.05921,-0.05263,
-0.04605,-0.03947,-0.03289,-0.02632,-0.01974,-0.01316,-0.00658,0,0.00658,
0.01316,0.01974,0.02632,0.03289,0.03947,0.04605,0.05263,0.05921,0.06579,
0.07237,0.07895,0.08553,0.09211,0.09868,0.10526,0.11184,0.11842,0.125,
0.13158,0.13816,0.14474,0.15132,0.15789,0.16447,0.17105,0.17763,0.18421,
0.19079,0.19737,0.20395,0.21053,0.21711,0.22368,0.23026,0.23684,0.24342,
0.25,0.25658,0.26316,0.26974,0.27632,0.28289,0.28947,0.29605,0.30263,
0.30921,0.31579,0.32237,0.32895,0.33553,0.34211,0.34868,0.35526,0.36184,
0.36842,0.375,0.38158,0.38816,0.39474,0.40132,0.40789,0.41447,0.42105,
0.42763,0.43421,0.44079,0.44737,0.45395,0.46053,0.46711,0.47368,0.48026,
0.48684,0.49342,0.5,0.50658,0.51316,0.51974,0.52632,0.53289,0.53947,
0.54605,0.55263,0.55921,0.56579,0.57237,0.57895,0.58553,0.59211,0.59868,
0.60526,0.61184,0.61842,0.625,0.63158,0.63816,0.64474,0.65132,0.65789,
0.66447,0.67105,0.67763,0.68421,0.69079,0.69737,0.70395,0.71053,0.71711,
0.72368,0.73026,0.73684,0.74342,0.75,0.75658,0.76316,0.76974,0.77632,
0.78289,0.78947,0.79605,0.80263,0.80921,0.81579,0.82237,0.82895,0.83553,
0.84211,0.84868,0.85526,0.86184,0.86842,0.875,0.88158,0.88816,0.89474,
0.90132,0.90789,0.91447,0.92105,0.92763,0.93421,0.94079,0.94737,0.95395,
0.96053,0.96711,0.97368,0.98026,0.98684,0.99342,1)

a2 <- c(1,0.98033,0.96079,0.94137,0.92209,0.90294,0.88392,0.86502,0.84626,0.82763,
0.80912,0.79075,0.77251,0.75439,0.73641,0.71856,0.70083,0.68324,0.66577,
0.64844,0.63123,0.61416,0.59721,0.5804,0.56371,0.54716,0.53073,0.51443,
0.49827,0.48223,0.46633,0.45055,0.4349,0.41939,0.404,0.38874,0.37361,
0.35862,0.34375,0.32901,0.3144,0.29993,0.28558,0.27136,0.25727,0.24331,
0.22948,0.21579,0.20222,0.18878,0.17547,0.16229,0.14924,0.13632,0.12353,
0.11087,0.09834,0.08594,0.07367,0.06153,0.04952,0.03763,0.02588,0.01426,
0.00277,-0.00859,-0.01982,-0.03093,-0.0419,-0.05274,-0.06345,-0.07403,
-0.08449,-0.09481,-0.105,-0.11507,-0.125,-0.1348,-0.14448,-0.15402,-0.16343,
-0.17272,-0.18187,-0.1909,-0.19979,-0.20856,-0.21719,-0.2257,-0.23407,
-0.24232,-0.25043,-0.25842,-0.26627,-0.274,-0.2816,-0.28906,-0.2964,
-0.30361,-0.31068,-0.31763,-0.32445,-0.33113,-0.33769,-0.34412,-0.35042,
-0.35658,-0.36262,-0.36853,-0.37431,-0.37996,-0.38547,-0.39086,-0.39612,
-0.40125,-0.40625,-0.41112,-0.41586,-0.42047,-0.42495,-0.4293,-0.43352,
-0.43761,-0.44157,-0.4454,-0.4491,-0.45267,-0.45611,-0.45942,-0.4626,
-0.46566,-0.46858,-0.47137,-0.47403,-0.47656,-0.47896,-0.48124,-0.48338,
-0.48539,-0.48727,-0.48903,-0.49065,-0.49214,-0.49351,-0.49474,-0.49584,
-0.49682,-0.49766,-0.49838,-0.49896,-0.49942,-0.49974,-0.49994,-0.5,
-0.49994,-0.49974,-0.49942,-0.49896,-0.49838,-0.49766,-0.49682,-0.49584,
-0.49474,-0.49351,-0.49214,-0.49065,-0.48903,-0.48727,-0.48539,-0.48338,
-0.48124,-0.47896,-0.47656,-0.47403,-0.47137,-0.46858,-0.46566,-0.4626,
-0.45942,-0.45611,-0.45267,-0.4491,-0.4454,-0.44157,-0.43761,-0.43352,
-0.4293,-0.42495,-0.42047,-0.41586,-0.41112,-0.40625,-0.40125,-0.39612,
-0.39086,-0.38547,-0.37996,-0.37431,-0.36853,-0.36262,-0.35658,-0.35042,
-0.34412,-0.33769,-0.33113,-0.32445,-0.31763,-0.31068,-0.30361,-0.2964,
-0.28906,-0.2816,-0.274,-0.26627,-0.25842,-0.25043,-0.24232,-0.23407,
-0.2257,-0.21719,-0.20856,-0.19979,-0.1909,-0.18187,-0.17272,-0.16343,
-0.15402,-0.14448,-0.1348,-0.125,-0.11507,-0.105,-0.09481,-0.08449,-0.07403,
-0.06345,-0.05274,-0.0419,-0.03093,-0.01982,-0.00859,0.00277,0.01426,
0.02588,0.03763,0.04952,0.06153,0.07367,0.08594,0.09834,0.11087,0.12353,
0.13632,0.14924,0.16229,0.17547,0.18878,0.20222,0.21579,0.22948,0.24331,
0.25727,0.27136,0.28558,0.29993,0.3144,0.32901,0.34375,0.35862,0.37361,
0.38874,0.404,0.41939,0.4349,0.45055,0.46633,0.48223,0.49827,0.51443,
0.53073,0.54716,0.56371,0.5804,0.59721,0.61416,0.63123,0.64844,0.66577,
0.68324,0.70083,0.71856,0.73641,0.75439,0.77251,0.79075,0.80912,0.82763,
0.84626,0.86502,0.88392,0.90294,0.92209,0.94137,0.96079,0.98033,1)

a3 <- c(-1,-0.96085,-0.92235,-0.88448,-0.84725,-0.81066,-0.77469,-0.73935,-0.70462,
-0.67051,-0.63701,-0.60412,-0.57183,-0.54014,-0.50904,-0.47853,-0.44861,
-0.41926,-0.3905,-0.3623,-0.33468,-0.30762,-0.28111,-0.25517,-0.22977,
-0.20492,-0.18061,-0.15685,-0.13361,-0.11091,-0.08873,-0.06707,-0.04593,
-0.0253,-0.00517,0.01444,0.03356,0.05218,0.07031,0.08796,0.10512,0.1218,
0.13801,0.15375,0.16902,0.18383,0.19819,0.21209,0.22554,0.23855,0.25112,
0.26326,0.27496,0.28623,0.29709,0.30752,0.31754,0.32715,0.33635,0.34515,
0.35356,0.36157,0.36919,0.37643,0.38329,0.38977,0.39589,0.40163,0.40701,
0.41203,0.4167,0.42102,0.42499,0.42862,0.43191,0.43487,0.4375,0.43981,
0.44179,0.44346,0.44482,0.44587,0.44661,0.44706,0.44721,0.44708,0.44665,
0.44595,0.44496,0.44371,0.44218,0.44039,0.43834,0.43603,0.43347,0.43066,
0.42761,0.42432,0.4208,0.41704,0.41306,0.40886,0.40444,0.3998,0.39496,
0.38991,0.38466,0.37921,0.37357,0.36774,0.36173,0.35554,0.34918,0.34264,
0.33594,0.32907,0.32205,0.31487,0.30755,0.30007,0.29246,0.28471,0.27683,
0.26882,0.26069,0.25244,0.24407,0.23559,0.227,0.21831,0.20953,0.20064,
0.19167,0.18262,0.17348,0.16427,0.15498,0.14562,0.1362,0.12673,0.11719,
0.10761,0.09797,0.0883,0.07858,0.06883,0.05906,0.04925,0.03943,0.02959,
0.01973,0.00987,0,-0.00987,-0.01973,-0.02959,-0.03943,-0.04925,-0.05906,
-0.06883,-0.07858,-0.0883,-0.09797,-0.10761,-0.11719,-0.12673,-0.1362,
-0.14562,-0.15498,-0.16427,-0.17348,-0.18262,-0.19167,-0.20064,-0.20953,
-0.21831,-0.227,-0.23559,-0.24407,-0.25244,-0.26069,-0.26882,-0.27683,
-0.28471,-0.29246,-0.30007,-0.30755,-0.31487,-0.32205,-0.32907,-0.33594,
-0.34264,-0.34918,-0.35554,-0.36173,-0.36774,-0.37357,-0.37921,-0.38466,
-0.38991,-0.39496,-0.3998,-0.40444,-0.40886,-0.41306,-0.41704,-0.4208,
-0.42432,-0.42761,-0.43066,-0.43347,-0.43603,-0.43834,-0.44039,-0.44218,
-0.44371,-0.44496,-0.44595,-0.44665,-0.44708,-0.44721,-0.44706,-0.44661,
-0.44587,-0.44482,-0.44346,-0.44179,-0.43981,-0.4375,-0.43487,-0.43191,
-0.42862,-0.42499,-0.42102,-0.4167,-0.41203,-0.40701,-0.40163,-0.39589,
-0.38977,-0.38329,-0.37643,-0.36919,-0.36157,-0.35356,-0.34515,-0.33635,
-0.32715,-0.31754,-0.30752,-0.29709,-0.28623,-0.27496,-0.26326,-0.25112,
-0.23855,-0.22554,-0.21209,-0.19819,-0.18383,-0.16902,-0.15375,-0.13801,
-0.1218,-0.10512,-0.08796,-0.07031,-0.05218,-0.03356,-0.01444,0.00517,
0.0253,0.04593,0.06707,0.08873,0.11091,0.13361,0.15685,0.18061,0.20492,
0.22977,0.25517,0.28111,0.30762,0.33468,0.3623,0.3905,0.41926,0.44861,
0.47853,0.50904,0.54014,0.57183,0.60412,0.63701,0.67051,0.70462,0.73935,
0.77469,0.81066,0.84725,0.88448,0.92235,0.96085,1)

a4 <- c(1,0.93518,0.87228,0.81126,0.75211,0.69478,0.63926,0.5855,0.53349,0.4832,
0.43459,0.38764,0.34232,0.2986,0.25646,0.21587,0.1768,0.13923,0.10312,
0.06845,0.0352,0.00334,-0.02716,-0.05632,-0.08417,-0.11074,-0.13604,
-0.1601,-0.18295,-0.20462,-0.22512,-0.24448,-0.26273,-0.27988,-0.29597,
-0.31101,-0.32503,-0.33805,-0.3501,-0.36119,-0.37135,-0.3806,-0.38896,
-0.39646,-0.40312,-0.40895,-0.41398,-0.41823,-0.42172,-0.42447,-0.4265,
-0.42784,-0.42849,-0.42849,-0.42784,-0.42658,-0.42472,-0.42227,-0.41926,
-0.41571,-0.41163,-0.40704,-0.40197,-0.39642,-0.39041,-0.38397,-0.37711,
-0.36985,-0.3622,-0.35418,-0.34581,-0.3371,-0.32807,-0.31874,-0.30911,
-0.29922,-0.28906,-0.27866,-0.26804,-0.25719,-0.24615,-0.23493,-0.22353,
-0.21198,-0.20028,-0.18845,-0.1765,-0.16445,-0.15231,-0.14009,-0.12781,
-0.11547,-0.10309,-0.09069,-0.07826,-0.06583,-0.0534,-0.04099,-0.0286,
-0.01626,-0.00396,0.00828,0.02045,0.03254,0.04455,0.05645,0.06825,0.07993,
0.09149,0.10291,0.11419,0.12532,0.13629,0.14709,0.15771,0.16816,0.17841,
0.18847,0.19832,0.20796,0.21739,0.22659,0.23556,0.24429,0.25279,0.26103,
0.26902,0.27676,0.28423,0.29143,0.29836,0.30502,0.31139,0.31747,0.32327,
0.32878,0.33399,0.3389,0.3435,0.3478,0.3518,0.35548,0.35885,0.36191,
0.36465,0.36707,0.36917,0.37095,0.37241,0.37354,0.37435,0.37484,0.375,
0.37484,0.37435,0.37354,0.37241,0.37095,0.36917,0.36707,0.36465,0.36191,
0.35885,0.35548,0.3518,0.3478,0.3435,0.3389,0.33399,0.32878,0.32327,
0.31747,0.31139,0.30502,0.29836,0.29143,0.28423,0.27676,0.26902,0.26103,
0.25279,0.24429,0.23556,0.22659,0.21739,0.20796,0.19832,0.18847,0.17841,
0.16816,0.15771,0.14709,0.13629,0.12532,0.11419,0.10291,0.09149,0.07993,
0.06825,0.05645,0.04455,0.03254,0.02045,0.00828,-0.00396,-0.01626,-0.0286,
-0.04099,-0.0534,-0.06583,-0.07826,-0.09069,-0.10309,-0.11547,-0.12781,
-0.14009,-0.15231,-0.16445,-0.1765,-0.18845,-0.20028,-0.21198,-0.22353,
-0.23493,-0.24615,-0.25719,-0.26804,-0.27866,-0.28906,-0.29922,-0.30911,
-0.31874,-0.32807,-0.3371,-0.34581,-0.35418,-0.3622,-0.36985,-0.37711,
-0.38397,-0.39041,-0.39642,-0.40197,-0.40704,-0.41163,-0.41571,-0.41926,
-0.42227,-0.42472,-0.42658,-0.42784,-0.42849,-0.42849,-0.42784,-0.4265,
-0.42447,-0.42172,-0.41823,-0.41398,-0.40895,-0.40312,-0.39646,-0.38896,
-0.3806,-0.37135,-0.36119,-0.3501,-0.33805,-0.32503,-0.31101,-0.29597,
-0.27988,-0.26273,-0.24448,-0.22512,-0.20462,-0.18295,-0.1601,-0.13604,
-0.11074,-0.08417,-0.05632,-0.02716,0.00334,0.0352,0.06845,0.10312,
0.13923,0.1768,0.21587,0.25646,0.2986,0.34232,0.38764,0.43459,0.4832,
0.53349,0.5855,0.63926,0.69478,0.75211,0.81126,0.87228,0.93518,1)

a5 <- c(-1,-0.90357,-0.81156,-0.72387,-0.64036,-0.56094,-0.48549,-0.41389,-0.34605,
-0.28185,-0.22119,-0.16396,-0.11007,-0.05941,-0.01188,0.0326,0.07414,
0.11283,0.14877,0.18203,0.21272,0.24092,0.26671,0.29018,0.31141,0.33048,
0.34747,0.36247,0.37554,0.38677,0.39622,0.40397,0.41009,0.41465,0.41772,
0.41936,0.41964,0.41863,0.41638,0.41296,0.40843,0.40285,0.39627,0.38875,
0.38035,0.37112,0.3611,0.35036,0.33895,0.3269,0.31427,0.30111,0.28746,
0.27336,0.25886,0.24399,0.22881,0.21334,0.19762,0.1817,0.16561,0.14939,
0.13306,0.11666,0.10022,0.08377,0.06735,0.05098,0.03468,0.0185,0.00244,
-0.01346,-0.02919,-0.04471,-0.06,-0.07506,-0.08984,-0.10435,-0.11855,
-0.13243,-0.14598,-0.15917,-0.172,-0.18444,-0.1965,-0.20814,-0.21937,
-0.23017,-0.24053,-0.25045,-0.2599,-0.2689,-0.27742,-0.28546,-0.29303,
-0.3001,-0.30668,-0.31276,-0.31835,-0.32343,-0.32801,-0.33209,-0.33566,
-0.33872,-0.34129,-0.34334,-0.3449,-0.34596,-0.34653,-0.3466,-0.34618,
-0.34528,-0.3439,-0.34204,-0.33972,-0.33694,-0.3337,-0.33001,-0.32589,
-0.32133,-0.31635,-0.31095,-0.30515,-0.29895,-0.29237,-0.28541,-0.27808,
-0.2704,-0.26238,-0.25403,-0.24535,-0.23637,-0.22709,-0.21753,-0.20769,
-0.1976,-0.18726,-0.1767,-0.16591,-0.15492,-0.14374,-0.13239,-0.12087,
-0.10921,-0.09741,-0.0855,-0.07348,-0.06137,-0.04918,-0.03694,-0.02465,
-0.01233,0,0.01233,0.02465,0.03694,0.04918,0.06137,0.07348,0.0855,0.09741,
0.10921,0.12087,0.13239,0.14374,0.15492,0.16591,0.1767,0.18726,0.1976,
0.20769,0.21753,0.22709,0.23637,0.24535,0.25403,0.26238,0.2704,0.27808,
0.28541,0.29237,0.29895,0.30515,0.31095,0.31635,0.32133,0.32589,0.33001,
0.3337,0.33694,0.33972,0.34204,0.3439,0.34528,0.34618,0.3466,0.34653,
0.34596,0.3449,0.34334,0.34129,0.33872,0.33566,0.33209,0.32801,0.32343,
0.31835,0.31276,0.30668,0.3001,0.29303,0.28546,0.27742,0.2689,0.2599,
0.25045,0.24053,0.23017,0.21937,0.20814,0.1965,0.18444,0.172,0.15917,
0.14598,0.13243,0.11855,0.10435,0.08984,0.07506,0.06,0.04471,0.02919,
0.01346,-0.00244,-0.0185,-0.03468,-0.05098,-0.06735,-0.08377,-0.10022,
-0.11666,-0.13306,-0.14939,-0.16561,-0.1817,-0.19762,-0.21334,-0.22881,
-0.24399,-0.25886,-0.27336,-0.28746,-0.30111,-0.31427,-0.3269,-0.33895,
-0.35036,-0.3611,-0.37112,-0.38035,-0.38875,-0.39627,-0.40285,-0.40843,
-0.41296,-0.41638,-0.41863,-0.41964,-0.41936,-0.41772,-0.41465,-0.41009,
-0.40397,-0.39622,-0.38677,-0.37554,-0.36247,-0.34747,-0.33048,-0.31141,
-0.29018,-0.26671,-0.24092,-0.21272,-0.18203,-0.14877,-0.11283,-0.07414,
-0.0326,0.01188,0.05941,0.11007,0.16396,0.22119,0.28185,0.34605,0.41389,
0.48549,0.56094,0.64036,0.72387,0.81156,0.90357,1)
}

if (std == "270"){

std <- 270

a0 <- rep(1, times = 270)

a1 <- c(-1,-0.99257,-0.98513,-0.9777,-0.97026,-0.96283,-0.95539,-0.94796,-0.94052,
-0.93309,-0.92565,-0.91822,-0.91078,-0.90335,-0.89591,-0.88848,-0.88104,
-0.87361,-0.86617,-0.85874,-0.8513,-0.84387,-0.83643,-0.829,-0.82156,
-0.81413,-0.80669,-0.79926,-0.79182,-0.78439,-0.77695,-0.76952,-0.76208,
-0.75465,-0.74721,-0.73978,-0.73234,-0.72491,-0.71747,-0.71004,-0.7026,
-0.69517,-0.68773,-0.6803,-0.67286,-0.66543,-0.65799,-0.65056,-0.64312,
-0.63569,-0.62825,-0.62082,-0.61338,-0.60595,-0.59851,-0.59108,-0.58364,
-0.57621,-0.56877,-0.56134,-0.5539,-0.54647,-0.53903,-0.5316,-0.52416,
-0.51673,-0.50929,-0.50186,-0.49442,-0.48699,-0.47955,-0.47212,-0.46468,
-0.45725,-0.44981,-0.44238,-0.43494,-0.42751,-0.42007,-0.41264,-0.4052,
-0.39777,-0.39033,-0.3829,-0.37546,-0.36803,-0.36059,-0.35316,-0.34572,
-0.33829,-0.33086,-0.32342,-0.31599,-0.30855,-0.30112,-0.29368,-0.28625,
-0.27881,-0.27138,-0.26394,-0.25651,-0.24907,-0.24164,-0.2342,-0.22677,
-0.21933,-0.2119,-0.20446,-0.19703,-0.18959,-0.18216,-0.17472,-0.16729,
-0.15985,-0.15242,-0.14498,-0.13755,-0.13011,-0.12268,-0.11524,-0.10781,
-0.10037,-0.09294,-0.0855,-0.07807,-0.07063,-0.0632,-0.05576,-0.04833,
-0.04089,-0.03346,-0.02602,-0.01859,-0.01115,-0.00372,0.00372,0.01115,
0.01859,0.02602,0.03346,0.04089,0.04833,0.05576,0.0632,0.07063,0.07807,
0.0855,0.09294,0.10037,0.10781,0.11524,0.12268,0.13011,0.13755,0.14498,
0.15242,0.15985,0.16729,0.17472,0.18216,0.18959,0.19703,0.20446,0.2119,
0.21933,0.22677,0.2342,0.24164,0.24907,0.25651,0.26394,0.27138,0.27881,
0.28625,0.29368,0.30112,0.30855,0.31599,0.32342,0.33086,0.33829,0.34572,
0.35316,0.36059,0.36803,0.37546,0.3829,0.39033,0.39777,0.4052,0.41264,
0.42007,0.42751,0.43494,0.44238,0.44981,0.45725,0.46468,0.47212,0.47955,
0.48699,0.49442,0.50186,0.50929,0.51673,0.52416,0.5316,0.53903,0.54647,
0.5539,0.56134,0.56877,0.57621,0.58364,0.59108,0.59851,0.60595,0.61338,
0.62082,0.62825,0.63569,0.64312,0.65056,0.65799,0.66543,0.67286,0.6803,
0.68773,0.69517,0.7026,0.71004,0.71747,0.72491,0.73234,0.73978,0.74721,
0.75465,0.76208,0.76952,0.77695,0.78439,0.79182,0.79926,0.80669,0.81413,
0.82156,0.829,0.83643,0.84387,0.8513,0.85874,0.86617,0.87361,0.88104,
0.88848,0.89591,0.90335,0.91078,0.91822,0.92565,0.93309,0.94052,0.94796,
0.95539,0.96283,0.97026,0.9777,0.98513,0.99257,1)

a2 <- c(1,0.97778,0.95572,0.93383,0.91211,0.89055,0.86916,0.84793,0.82687,0.80597,
0.78524,0.76468,0.74428,0.72405,0.70398,0.68408,0.66435,0.64478,0.62538,
0.60614,0.58707,0.56817,0.54943,0.53085,0.51244,0.4942,0.47613,0.45822,
0.44047,0.42289,0.40548,0.38823,0.37115,0.35424,0.33749,0.3209,0.30449,
0.28824,0.27215,0.25623,0.24047,0.22489,0.20946,0.19421,0.17912,0.16419,
0.14943,0.13484,0.12041,0.10615,0.09205,0.07812,0.06436,0.05076,0.03733,
0.02406,0.01096,-0.00198,-0.01475,-0.02735,-0.03979,-0.05206,-0.06416,
-0.0761,-0.08788,-0.09949,-0.11093,-0.12221,-0.13332,-0.14426,-0.15504,
-0.16566,-0.1761,-0.18638,-0.1965,-0.20645,-0.21624,-0.22585,-0.23531,
-0.24459,-0.25371,-0.26267,-0.27146,-0.28008,-0.28854,-0.29683,-0.30496,
-0.31292,-0.32071,-0.32834,-0.3358,-0.3431,-0.35023,-0.3572,-0.36399,
-0.37063,-0.3771,-0.3834,-0.38953,-0.3955,-0.40131,-0.40695,-0.41242,
-0.41773,-0.42287,-0.42784,-0.43265,-0.43729,-0.44177,-0.44608,-0.45023,
-0.45421,-0.45802,-0.46167,-0.46515,-0.46847,-0.47162,-0.47461,-0.47743,
-0.48008,-0.48257,-0.48489,-0.48704,-0.48903,-0.49086,-0.49252,-0.49401,
-0.49534,-0.4965,-0.49749,-0.49832,-0.49898,-0.49948,-0.49981,-0.49998,
-0.49998,-0.49981,-0.49948,-0.49898,-0.49832,-0.49749,-0.4965,-0.49534,
-0.49401,-0.49252,-0.49086,-0.48903,-0.48704,-0.48489,-0.48257,-0.48008,
-0.47743,-0.47461,-0.47162,-0.46847,-0.46515,-0.46167,-0.45802,-0.45421,
-0.45023,-0.44608,-0.44177,-0.43729,-0.43265,-0.42784,-0.42287,-0.41773,
-0.41242,-0.40695,-0.40131,-0.3955,-0.38953,-0.3834,-0.3771,-0.37063,
-0.36399,-0.3572,-0.35023,-0.3431,-0.3358,-0.32834,-0.32071,-0.31292,
-0.30496,-0.29683,-0.28854,-0.28008,-0.27146,-0.26267,-0.25371,-0.24459,
-0.23531,-0.22585,-0.21624,-0.20645,-0.1965,-0.18638,-0.1761,-0.16566,
-0.15504,-0.14426,-0.13332,-0.12221,-0.11093,-0.09949,-0.08788,-0.0761,
-0.06416,-0.05206,-0.03979,-0.02735,-0.01475,-0.00198,0.01096,0.02406,
0.03733,0.05076,0.06436,0.07812,0.09205,0.10615,0.12041,0.13484,0.14943,
0.16419,0.17912,0.19421,0.20946,0.22489,0.24047,0.25623,0.27215,0.28824,
0.30449,0.3209,0.33749,0.35424,0.37115,0.38823,0.40548,0.42289,0.44047,
0.45822,0.47613,0.4942,0.51244,0.53085,0.54943,0.56817,0.58707,0.60614,
0.62538,0.64478,0.66435,0.68408,0.70398,0.72405,0.74428,0.76468,0.78524,
0.80597,0.82687,0.84793,0.86916,0.89055,0.91211,0.93383,0.95572,0.97778,
1)

a3 <- c(-1,-0.9558,-0.91243,-0.86987,-0.82813,-0.78719,-0.74705,-0.70769,-0.66913,
-0.63135,-0.59433,-0.55809,-0.52261,-0.48788,-0.4539,-0.42067,-0.38817,
-0.3564,-0.32536,-0.29504,-0.26542,-0.23651,-0.20831,-0.18079,-0.15397,
-0.12782,-0.10235,-0.07755,-0.05341,-0.02993,-0.0071,0.01509,0.03664,
0.05756,0.07785,0.09752,0.11658,0.13503,0.15288,0.17014,0.1868,0.20289,
0.2184,0.23333,0.24771,0.26152,0.27479,0.28751,0.29968,0.31133,0.32245,
0.33305,0.34313,0.3527,0.36177,0.37035,0.37844,0.38604,0.39316,0.39981,
0.406,0.41173,0.417,0.42183,0.42621,0.43017,0.43369,0.43679,0.43947,
0.44175,0.44362,0.44509,0.44618,0.44687,0.44719,0.44714,0.44671,0.44593,
0.44479,0.44331,0.44148,0.43932,0.43682,0.43401,0.43087,0.42742,0.42367,
0.41962,0.41528,0.41065,0.40574,0.40056,0.3951,0.38939,0.38342,0.3772,
0.37073,0.36403,0.3571,0.34994,0.34257,0.33498,0.32718,0.31919,0.311,
0.30262,0.29406,0.28532,0.27642,0.26735,0.25812,0.24875,0.23923,0.22957,
0.21977,0.20985,0.19981,0.18966,0.1794,0.16904,0.15858,0.14803,0.1374,
0.12669,0.11591,0.10507,0.09416,0.08321,0.07221,0.06117,0.05009,0.03899,
0.02786,0.01673,0.00558,-0.00558,-0.01673,-0.02786,-0.03899,-0.05009,
-0.06117,-0.07221,-0.08321,-0.09416,-0.10507,-0.11591,-0.12669,-0.1374,
-0.14803,-0.15858,-0.16904,-0.1794,-0.18966,-0.19981,-0.20985,-0.21977,
-0.22957,-0.23923,-0.24875,-0.25812,-0.26735,-0.27642,-0.28532,-0.29406,
-0.30262,-0.311,-0.31919,-0.32718,-0.33498,-0.34257,-0.34994,-0.3571,
-0.36403,-0.37073,-0.3772,-0.38342,-0.38939,-0.3951,-0.40056,-0.40574,
-0.41065,-0.41528,-0.41962,-0.42367,-0.42742,-0.43087,-0.43401,-0.43682,
-0.43932,-0.44148,-0.44331,-0.44479,-0.44593,-0.44671,-0.44714,-0.44719,
-0.44687,-0.44618,-0.44509,-0.44362,-0.44175,-0.43947,-0.43679,-0.43369,
-0.43017,-0.42621,-0.42183,-0.417,-0.41173,-0.406,-0.39981,-0.39316,
-0.38604,-0.37844,-0.37035,-0.36177,-0.3527,-0.34313,-0.33305,-0.32245,
-0.31133,-0.29968,-0.28751,-0.27479,-0.26152,-0.24771,-0.23333,-0.2184,
-0.20289,-0.1868,-0.17014,-0.15288,-0.13503,-0.11658,-0.09752,-0.07785,
-0.05756,-0.03664,-0.01509,0.0071,0.02993,0.05341,0.07755,0.10235,0.12782,
0.15397,0.18079,0.20831,0.23651,0.26542,0.29504,0.32536,0.3564,0.38817,
0.42067,0.4539,0.48788,0.52261,0.55809,0.59433,0.63135,0.66913,0.70769,
0.74705,0.78719,0.82813,0.86987,0.91243,0.9558,1)

a4 <- c(1,0.92689,0.85622,0.78795,0.72205,0.65846,0.59714,0.53806,0.48118,0.42644,
0.37382,0.32327,0.27476,0.22823,0.18366,0.14101,0.10023,0.06129,0.02415,
-0.01123,-0.04488,-0.07685,-0.10716,-0.13585,-0.16297,-0.18854,-0.2126,
-0.23519,-0.25634,-0.27609,-0.29446,-0.3115,-0.32723,-0.34169,-0.35491,
-0.36693,-0.37777,-0.38748,-0.39607,-0.40358,-0.41004,-0.41549,-0.41995,
-0.42344,-0.42602,-0.42769,-0.42849,-0.42845,-0.42759,-0.42595,-0.42355,
-0.42042,-0.41659,-0.41208,-0.40692,-0.40113,-0.39474,-0.38778,-0.38027,
-0.37224,-0.36371,-0.3547,-0.34524,-0.33535,-0.32505,-0.31437,-0.30333,
-0.29196,-0.28026,-0.26828,-0.25601,-0.2435,-0.23075,-0.21779,-0.20464,
-0.19132,-0.17784,-0.16423,-0.1505,-0.13668,-0.12277,-0.1088,-0.09479,
-0.08075,-0.0667,-0.05266,-0.03864,-0.02465,-0.01072,0.00315,0.01693,
0.03062,0.04419,0.05764,0.07095,0.08411,0.09711,0.10993,0.12256,0.13499,
0.14721,0.1592,0.17096,0.18247,0.19373,0.20473,0.21545,0.22588,0.23602,
0.24586,0.25539,0.2646,0.27348,0.28203,0.29025,0.29811,0.30562,0.31277,
0.31956,0.32597,0.33201,0.33766,0.34294,0.34782,0.35231,0.3564,0.36009,
0.36338,0.36627,0.36874,0.37081,0.37246,0.3737,0.37453,0.37495,0.37495,
0.37453,0.3737,0.37246,0.37081,0.36874,0.36627,0.36338,0.36009,0.3564,
0.35231,0.34782,0.34294,0.33766,0.33201,0.32597,0.31956,0.31277,0.30562,
0.29811,0.29025,0.28203,0.27348,0.2646,0.25539,0.24586,0.23602,0.22588,
0.21545,0.20473,0.19373,0.18247,0.17096,0.1592,0.14721,0.13499,0.12256,
0.10993,0.09711,0.08411,0.07095,0.05764,0.04419,0.03062,0.01693,0.00315,
-0.01072,-0.02465,-0.03864,-0.05266,-0.0667,-0.08075,-0.09479,-0.1088,
-0.12277,-0.13668,-0.1505,-0.16423,-0.17784,-0.19132,-0.20464,-0.21779,
-0.23075,-0.2435,-0.25601,-0.26828,-0.28026,-0.29196,-0.30333,-0.31437,
-0.32505,-0.33535,-0.34524,-0.3547,-0.36371,-0.37224,-0.38027,-0.38778,
-0.39474,-0.40113,-0.40692,-0.41208,-0.41659,-0.42042,-0.42355,-0.42595,
-0.42759,-0.42845,-0.42849,-0.42769,-0.42602,-0.42344,-0.41995,-0.41549,
-0.41004,-0.40358,-0.39607,-0.38748,-0.37777,-0.36693,-0.35491,-0.34169,
-0.32723,-0.3115,-0.29446,-0.27609,-0.25634,-0.23519,-0.2126,-0.18854,
-0.16297,-0.13585,-0.10716,-0.07685,-0.04488,-0.01123,0.02415,0.06129,
0.10023,0.14101,0.18366,0.22823,0.27476,0.32327,0.37382,0.42644,0.48118,
0.53806,0.59714,0.65846,0.72205,0.78795,0.85622,0.92689,1)

a5 <- c(-1,-0.89135,-0.78833,-0.69078,-0.59853,-0.51141,-0.42927,-0.35195,-0.2793,
-0.21116,-0.14739,-0.08783,-0.03235,0.01919,0.06694,0.11103,0.15159,
0.18875,0.22264,0.25339,0.28111,0.30594,0.32798,0.34736,0.36418,0.37855,
0.39059,0.4004,0.40809,0.41375,0.41748,0.41939,0.41956,0.41809,0.41507,
0.41059,0.40472,0.39757,0.3892,0.37969,0.36913,0.35759,0.34514,0.33186,
0.3178,0.30305,0.28766,0.27171,0.25524,0.23833,0.22102,0.20337,0.18545,
0.16729,0.14896,0.1305,0.11195,0.09337,0.07479,0.05627,0.03783,0.01952,
0.00137,-0.01658,-0.03429,-0.05173,-0.06888,-0.08569,-0.10216,-0.11823,
-0.13391,-0.14915,-0.16393,-0.17824,-0.19206,-0.20537,-0.21814,-0.23037,
-0.24204,-0.25313,-0.26364,-0.27355,-0.28286,-0.29155,-0.29961,-0.30705,
-0.31386,-0.32003,-0.32555,-0.33044,-0.33467,-0.33827,-0.34122,-0.34352,
-0.34519,-0.34622,-0.34662,-0.3464,-0.34555,-0.34409,-0.34202,-0.33936,
-0.3361,-0.33227,-0.32787,-0.32292,-0.31742,-0.31139,-0.30484,-0.29778,
-0.29024,-0.28221,-0.27373,-0.2648,-0.25545,-0.24568,-0.23552,-0.22498,
-0.21408,-0.20285,-0.19129,-0.17943,-0.16729,-0.15488,-0.14224,-0.12937,
-0.11629,-0.10304,-0.08963,-0.07608,-0.0624,-0.04864,-0.0348,-0.0209,
-0.00697,0.00697,0.0209,0.0348,0.04864,0.0624,0.07608,0.08963,0.10304,
0.11629,0.12937,0.14224,0.15488,0.16729,0.17943,0.19129,0.20285,0.21408,
0.22498,0.23552,0.24568,0.25545,0.2648,0.27373,0.28221,0.29024,0.29778,
0.30484,0.31139,0.31742,0.32292,0.32787,0.33227,0.3361,0.33936,0.34202,
0.34409,0.34555,0.3464,0.34662,0.34622,0.34519,0.34352,0.34122,0.33827,
0.33467,0.33044,0.32555,0.32003,0.31386,0.30705,0.29961,0.29155,0.28286,
0.27355,0.26364,0.25313,0.24204,0.23037,0.21814,0.20537,0.19206,0.17824,
0.16393,0.14915,0.13391,0.11823,0.10216,0.08569,0.06888,0.05173,0.03429,
0.01658,-0.00137,-0.01952,-0.03783,-0.05627,-0.07479,-0.09337,-0.11195,
-0.1305,-0.14896,-0.16729,-0.18545,-0.20337,-0.22102,-0.23833,-0.25524,
-0.27171,-0.28766,-0.30305,-0.3178,-0.33186,-0.34514,-0.35759,-0.36913,
-0.37969,-0.3892,-0.39757,-0.40472,-0.41059,-0.41507,-0.41809,-0.41956,
-0.41939,-0.41748,-0.41375,-0.40809,-0.4004,-0.39059,-0.37855,-0.36418,
-0.34736,-0.32798,-0.30594,-0.28111,-0.25339,-0.22264,-0.18875,-0.15159,
-0.11103,-0.06694,-0.01919,0.03235,0.08783,0.14739,0.21116,0.2793,0.35195,
0.42927,0.51141,0.59853,0.69078,0.78833,0.89135,1)
}

# Create list of cows:
cows <- unique(dat[,1])

rval <- as.data.frame(matrix(rep(0, times = 9), nrow = 1))
names(rval) <- c("ckey", "lact", "llen", "vltd", "fltd", "pltd", "vstd", "fstd", "pstd")

# Convert fat and protein yields to kilograms (1 litre milk = 0.970264 kg):
dat$fat <- dat$fat * (dat$vol * 0.970264) 
dat$pro <- dat$pro * (dat$vol * 0.970264)

# Loop through each cows's records:
for(i in 1:length(cows)){

 # Select herd test records for this cow:
 id <- dat[,1] == cows[i]
 dat.tmp <- dat[id, 1:8]
 
 # Take each lactation in turn:
 lacts <- unique(dat.tmp[,3])
 
 for(j in 1:length(lacts)){
  id <- dat.tmp[,3] == lacts[j]
  dat.ump <- dat.tmp[id,1:8]

  # How many herd tests were there?
  ntest <- length(dat.ump[,4])

  # If there less than 4 herd test events, don't calculate anything:
  if (ntest < 4){
   vstd <- 0; fstd <- 0; pstd <- 0
   vltd <- 0; fltd <- 0; pltd <- 0
  }

  if (ntest >= 4){
  # Extract appropriate values of x on the basis of herd test days in milk:
  x0 <- a0[dat.ump[,4]]
  x1 <- a1[dat.ump[,4]]
  x2 <- a2[dat.ump[,4]]
  x3 <- a3[dat.ump[,4]]
  xmat <- cbind(x0, x1, x2, x3)
  txmat <- t(xmat)
  Xx <- txmat %*% xmat 
  
  Xy.vol <- txmat %*% dat.ump[,6]
  Xy.fat <- txmat %*% dat.ump[,7]
  Xy.pro <- txmat %*% dat.ump[,8]
 
  # Regression coefficients: 
  a.vol <- solve(Xx, Xy.vol)
  a.fat <- solve(Xx, Xy.fat)
  a.pro <- solve(Xx, Xy.pro)

  # Variable "flag" equals TRUE if there is no dry off date and FALSE otherwise:
  flag <- is.na(unique(dat.ump[,5]))

  # Lactation length and days in milk at last herd test:
  llen <- unique(dat.ump[,5])
  last.ht <- max(dat.ump[,4])

  # If dry (i.e. flag == FALSE) and llen <= std, calculate yield to dry off date:
  if (flag == FALSE & llen <= std){
   vstd <- sum((a0[1:std] * a.vol[1]) +  (a1[1:std] * a.vol[2]) + (a2[1:std] * a.vol[3]) + (a4[1:std] * a.vol[3]))
   fstd <- sum((a0[1:std] * a.fat[1]) +  (a1[1:std] * a.fat[2]) + (a2[1:std] * a.fat[3]) + (a4[1:std] * a.fat[3]))
   pstd <- sum((a0[1:std] * a.pro[1]) +  (a1[1:std] * a.pro[2]) + (a2[1:std] * a.pro[3]) + (a4[1:std] * a.pro[3]))

   vltd <- sum((a0[1:llen] * a.vol[1]) +  (a1[1:llen] * a.vol[2]) + (a2[1:llen] * a.vol[3]) + (a4[1:llen] * a.vol[3]))
   fltd <- sum((a0[1:llen] * a.fat[1]) +  (a1[1:llen] * a.fat[2]) + (a2[1:llen] * a.fat[3]) + (a4[1:llen] * a.fat[3]))
   pltd <- sum((a0[1:llen] * a.pro[1]) +  (a1[1:llen] * a.pro[2]) + (a2[1:llen] * a.pro[3]) + (a4[1:llen] * a.pro[3]))
   }
  
  # If dry (i.e. flag == FALSE) and llen > std, calculate yield to dry off date:
  if (flag == FALSE & llen > std){
   vstd <- sum((a0[1:std] * a.vol[1]) +  (a1[1:std] * a.vol[2]) + (a2[1:std] * a.vol[3]) + (a4[1:std] * a.vol[3]))
   fstd <- sum((a0[1:std] * a.fat[1]) +  (a1[1:std] * a.fat[2]) + (a2[1:std] * a.fat[3]) + (a4[1:std] * a.fat[3]))
   pstd <- sum((a0[1:std] * a.pro[1]) +  (a1[1:std] * a.pro[2]) + (a2[1:std] * a.pro[3]) + (a4[1:std] * a.pro[3]))

   d.xs <- llen - std
   v.xs <- (sum((a0[std]*a.vol[1])+(a1[std]*a.vol[2]) + (a2[std]*a.vol[3])+(a4[std]*a.vol[3]))) + (dat.ump[,6][dat.ump[,4]==last.ht]/2)
   f.xs <- (sum((a0[std]*a.fat[1])+(a1[std]*a.fat[2]) + (a2[std]*a.fat[3])+(a4[std]*a.fat[3]))) + (dat.ump[,7][dat.ump[,4]==last.ht]/2)
   p.xs <- (sum((a0[std]*a.pro[1])+(a1[std]*a.pro[2]) + (a2[std]*a.pro[3])+(a4[std]*a.pro[3]))) + (dat.ump[,8][dat.ump[,4]==last.ht]/2)
   
   vltd <- vstd + (v.xs * d.xs)
   fltd <- fstd + (f.xs * d.xs)
   pltd <- pstd + (p.xs * d.xs)
   }
  
  # If lactating (i.e. flag == TRUE) and dim at last herd test <= std, calculate yield to last herd test:
  if (flag == TRUE & last.ht <= std){
   vstd <- sum((a0[1:std] * a.vol[1]) +  (a1[1:std] * a.vol[2]) + (a2[1:std] * a.vol[3]) + (a4[1:std] * a.vol[3]))
   fstd <- sum((a0[1:std] * a.fat[1]) +  (a1[1:std] * a.fat[2]) + (a2[1:std] * a.fat[3]) + (a4[1:std] * a.fat[3]))
   pstd <- sum((a0[1:std] * a.pro[1]) +  (a1[1:std] * a.pro[2]) + (a2[1:std] * a.pro[3]) + (a4[1:std] * a.pro[3]))

   vltd <- sum((a0[1:last.ht] * a.vol[1]) + (a1[1:last.ht] * a.vol[2]) + (a2[1:last.ht] * a.vol[3]) + (a4[1:last.ht] * a.vol[3]))
   fltd <- sum((a0[1:last.ht] * a.fat[1]) + (a1[1:last.ht] * a.fat[2]) + (a2[1:last.ht] * a.fat[3]) + (a4[1:last.ht] * a.fat[3]))
   pltd <- sum((a0[1:last.ht] * a.pro[1]) + (a1[1:last.ht] * a.pro[2]) + (a2[1:last.ht] * a.pro[3]) + (a4[1:last.ht] * a.pro[3]))
   }

  # If lactating (i.e. flag == TRUE) and dim at last herd test > std, calculate yield to last herd test:
  if (flag == TRUE & last.ht > std){
   vstd <- sum((a0[1:std] * a.vol[1]) +  (a1[1:std] * a.vol[2]) + (a2[1:std] * a.vol[3]) + (a4[1:std] * a.vol[3]))
   fstd <- sum((a0[1:std] * a.fat[1]) +  (a1[1:std] * a.fat[2]) + (a2[1:std] * a.fat[3]) + (a4[1:std] * a.fat[3]))
   pstd <- sum((a0[1:std] * a.pro[1]) +  (a1[1:std] * a.pro[2]) + (a2[1:std] * a.pro[3]) + (a4[1:std] * a.pro[3]))

   d.xs <- last.ht - std
   v.xs <- (sum((a0[std]*a.vol[1])+(a1[std]*a.vol[2]) + (a2[std]*a.vol[3])+(a4[std]*a.vol[3]))) + (dat.ump[,6][dat.ump[,4]==last.ht]/2)
   f.xs <- (sum((a0[std]*a.fat[1])+(a1[std]*a.fat[2]) + (a2[std]*a.fat[3])+(a4[std]*a.fat[3]))) + (dat.ump[,7][dat.ump[,4]==last.ht]/2)
   p.xs <- (sum((a0[std]*a.pro[1])+(a1[std]*a.pro[2]) + (a2[std]*a.pro[3])+(a4[std]*a.pro[3]))) + (dat.ump[,8][dat.ump[,4]==last.ht]/2)
   
   vltd <- vstd + (v.xs * d.xs)
   fltd <- fstd + (f.xs * d.xs)
   pltd <- pstd + (p.xs * d.xs)
   }
  }
  ckey <- unique(dat.ump[,1])
  lact <- unique(dat.ump[,3])
  rval.tmp <- round(as.data.frame(cbind(ckey, lact, llen, vltd, fltd, pltd, vstd, fstd, pstd)), digits = 0)
  rval <- rbind(rval, rval.tmp)
 }
}
 rval <- as.data.frame(rval[-1,], row.names = NULL)
 print(rval)
}