tests/testthat/test-loglikelihood-covariates.R

context("Contrast and log-likelihood with covariates")


test_that("Weibull+ARAInf+2cov",{
			 DataC<-data.frame(System=c(rep(1,5),rep(2,4),rep(3,1),rep(4,5)),Time=c(0.800,2.646,3.190,3.916,4.109,0.910,1.127,1.245, 1.349, 0.541,1.397,1.463,2.406,2.506,3.159),Type=rep(-1,15))
			 Cov<-data.frame(cov1=c(4.336,5.615,4.770,4.655),cov2=c(0,0,0,1))
			 mle<-mle.vam(System&Time&Type~(ARAInf(0.8)|Weibull(0.15,2.3|0.6*cov1-0.9*cov2)),data=DataC,data.covariates = Cov)
			 theta<-c(0.15,2.3,0.8,0.6,-0.9)
 lnL=-10.611371878394
 dlnL=c(-52.9524008304927,-7.48846756301192,-5.00289293789225,-33.9797994354049,1.30146762357199)
 d2lnL=matrix(c(-666.666666666667,-39.4713839842167,174.067812605015,-707.765329569366,-24.6568825095201,-39.4713839842167,-8.55601730656781,-5.06266541120669,-24.5273121946206,-0.526160464419239,174.067812605015,-5.06266541120669,-80.3899677195699,117.185366639259,2.57743265786856,-707.765329569366,-24.5273121946206,117.185366639259,-496.497440827701,-17.2166682122724,-24.6568825095201,-0.526160464419239,2.57743265786856,-17.2166682122724,-3.69853237642801),nrow=5,byrow=TRUE)
 C=-9.04286047358826
 dC=c(-5.43870776710372,-14.0422821871295,2.77464824418773,2.58190629480419)
 d2C=matrix(c(-5.57662804602105,-16.0009942136532,1.87634470340362,0.280018217910437,-16.0009942136532,-41.6377558310791,-2.37702341629769,-1.06679604415864,1.87634470340362,-2.37702341629769,-3.42265694392171,-0.0668438872760557,0.280018217910437,-1.06679604415864,-0.0668438872760557,-2.02828189405531),nrow=4,byrow=TRUE)
 expect_that(logLik.mle.vam(mle) ,equals(lnL,tolerance=0.00000000000001))
	expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
  expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
  expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
  expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
  expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
  expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
}
)

test_that("Weibull+ARA1+cov+1syst",{
				DataC<-data.frame(System=c(rep(1,5)),Time=c(0.800,2.646,3.190,3.916,4.109),Type=rep(-1,5))
			 Cov<-data.frame(cov1=c(4.336))
			 mle<-mle.vam(System&Time&Type~(ARA1(0.8)|Weibull(0.15,2.3|0.6*cov1)),data=DataC,data.covariates = Cov)
			 theta<-c(0.15,2.3,0.8,0.6)
			 lnL=-8.52877195097165
dlnL=c(-81.9645913694061,-6.07046578028585,23.3314139214845,-53.3097702266617)
d2lnL=matrix(c(-222.222222222222,-60.495710436341,235.401846082455,-499.931801511078,-60.495710436341,-5.77322570679121,30.9819321822151,-39.3464100677962,235.401846082455,30.9819321822151,-66.6000012977875,153.105360692028,-499.931801511078,-39.3464100677962,153.105360692028,-325.155643702805),nrow=4,byrow=TRUE)
C=-2.43889089898127
dC=c(0.38043894974369,-1.77044654759699,-3.5527136788005e-15)
d2C=matrix(c(-0.964497262058924,-2.94973113824478,0,-2.94973113824478,-24.1032000234141,0,0,0,0),nrow=3,byrow=TRUE)
expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
}
)

test_that("loglinear+GQRARA3+1cov",{
DataC<-data.frame(System=c(rep(1,5),rep(2,4),rep(3,1),rep(4,5)),Time=c(0.800,2.646,3.190,3.916,4.109,0.910,1.127,1.245, 1.349, 0.541,1.397,1.463,2.406,2.506,3.159),Type=rep(-1,15))
 Cov<-data.frame(cov1=c(4.336,5.615,4.770,4.655))
 mle<-mle.vam(System&Time&Type~(GQR_ARAm(0.9,0.6|log,3)|LogLinear(0.01,0.8|0.6*cov1)),data=DataC,data.covariates = Cov)
 theta<-c(10,0.8,0.9,0.6,0.05)
 lnL=-156.309169533704
 dlnL=c(-18.8469478697017,-165.776432808115,-201.625671279627,123.213988284382,-865.979294381439)
 d2lnL=matrix(c(-0.15,-17.8721369503378,-22.2292998245236,13.581946596117,-93.8164294381439,-17.8721369503378,-206.482527176374,-318.907779510921,284.196589618172,-805.979097834773,-22.2292998245236,-318.907779510921,-254.20001606188,296.192921955,-993.304604477443,13.581946596117,284.196589618172,296.192921955,-345.768395975088,605.804302391538,-93.8164294381439,-805.979097834773,-993.304604477443,605.804302391538,-4358.96682220658),nrow=5,byrow=TRUE)
 C=-6.9516779971749
 dC=c(-0.230604716245915,4.27963582376563,-2.59271309674429,3.02246758560939)
 d2C=matrix(c(-3.64913257451202,-1.59292232002841,-2.43884692315837,1.33256371243673,-1.59292232002841,-14.6491411100637,5.3104285778739,2.3334095733029,-2.43884692315837,5.3104285778739,3.28781023351844,-1.50663226315177,1.33256371243673,2.3334095733029,-1.50663226315177,-2.4509027895092),nrow=4,byrow=TRUE)
 expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
 expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
}
)

test_that("Weibull+GQRARA3+PM_AGAN+2cov",{
			DataC<-data.frame(System=c(rep(1,5),rep(2,4),rep(3,1),rep(4,5)),Time=c(0.800,2.646,3.190,3.916,4.109,0.910,1.127,1.245, 1.349, 0.541,1.397,1.463,2.406,2.506,3.159),Type=c(1,-1,1,1,0, -1,-1,1,-1, 0, 1,-1,-1,1,1))
			 Cov<-data.frame(cov1=c(4.336,5.615,4.770,4.655),cov2=c(1,0,0,1))
			 mle<-mle.vam(System&Time&Type~(GQR_ARAm(0.9,0.6|log,3)|Weibull(0.1,2.1|0.6*cov1+0.7*cov2))&(AGAN()),data=DataC,data.covariates = Cov)
			 theta<-c(0.1,2.5,0.9,0.3,0.6,-0.9)
			 lnL=-10.7825703110304
 dlnL=c(-48.3485403654731,-4.81807237872075,-1.06819299760573,3.99004211321699,-22.0525403274032,-3.83179826270904)
 d2lnL=matrix(c(-600,-24.9907943270565,-39.8377159829799,57.2481109937894,-525.435403274032,-68.3179826270904,-24.9907943270565,-2.95404657311797,-1.0146194427903,4.7267032956704,-10.6476306501269,-2.89212755325284,-39.8377159829799,-1.0146194427903,-4.40339100506147,6.26245971708459,-19.5765879978775,-2.48100912650979,57.2481109937894,4.7267032956704,6.26245971708459,-7.64915935008439,28.8333953267537,2.6575443253662,-525.435403274032,-10.6476306501269,-19.5765879978775,28.8333953267537,-258.198848960839,-30.3846558030683,-68.3179826270904,-2.89212755325284,-2.48100912650979,2.6575443253662,-30.3846558030683,-6.83179826270904),nrow=6,byrow=TRUE)
 C=-9.4937684396386
 dC=c(-3.70290421474614,0.709491697337453,1.43545061074362,1.39404669532198,-0.783234129226588)
 d2C=matrix(c(-1.74503828677723,0.356094324870721,1.37020944094029,0.814945326219616,-0.728958848122487,0.356094324870721,-2.98149970164565,2.47186711348157,-0.142494934392076,0.0171193736473143,1.37020944094029,2.47186711348157,-3.35219640541417,0.593068508548029,-0.527283040738493,0.814945326219616,-0.142494934392076,0.593068508548029,-1.87692001456765,1.52070163225436,-0.728958848122487,0.0171193736473143,-0.527283040738493,1.52070163225436,-1.39775738313578),nrow=5,byrow=TRUE)
 expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
 expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
})

test_that("Weibull+ABAO+PM_AGAN+2cov",{
			DataC<-data.frame(System=c(rep(1,5),rep(2,4),rep(3,1),rep(4,5)),Time=c(0.800,2.646,3.190,3.916,4.109,0.910,1.127,1.245, 1.349, 0.541,1.397,1.463,2.406,2.506,3.159),Type=c(1,-1,1,1,0, -1,-1,1,-1, 0, 1,-1,-1,1,1))
         Cov<-data.frame(cov1=c(4.336,5.615,4.770,4.655),cov2=c(1,0,0,1))
         mle<-mle.vam(System&Time&Type~(ABAO()|Weibull(0.1,2.1|0.6*cov1+0.7*cov2))&(AGAN()),data=DataC,data.covariates = Cov)
         theta<-c(0.1,2.5,0.6,-0.9)
				 lnL=-12.6662572894128
 dlnL=c(-74.4173789705649,-7.33918738873799,-35.2305026092892,-5.03081291162338)
 d2lnL=matrix(c(-600,-54.0495190498004,-657.215026092892,-80.3081291162337,-54.0495190498004,-5.3579333086138,-24.8748617948241,-4.5583115895231,-657.215026092892,-24.8748617948241,-325.792993578815,-35.6574324618493,-80.3081291162337,-4.5583115895231,-35.6574324618493,-8.03081291162338),nrow=4,byrow=TRUE)
 C=-10.064150384301
 dC=c(-4.34684847162031,1.15483687320025,-0.584720803125606)
 d2C=matrix(c(-1.95299233396594,0.692721438207464,-0.593273434174727,0.692721438207464,-1.98941834703072,1.61055385924745,-0.593273434174727,1.61055385924745,-1.44301693039869),nrow=3,byrow=TRUE)
 expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
 expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
})

test_that("Weibull+GQR_ARAm+PM_ARA1+PM_GQR_ARAm+2cov",{
			DataC<-data.frame(System=c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4),Time=c(0.79,1.583,1.761,2,2.652,2.705,3.445,4,4.57,4.901,0.069,0.477,0.682,1.004,1.327,1.704,1.83,1.848,0.735,1.004,1.408,1.874,2,2.466,2.514,2.919,3.376,1.02,1.401,2,2.216,2.752,3.426,3.592,4,4.35,4.82),Type=c(-1,-1,-1,2,1,-1,1,2,-1,-1,-1,-1,-1,-1,1,-1,-1,-1,1,-1,-1,1,2,-1,-1,1,1,-1,-1,2,-1,1,1,-1,2,-1,1))
			 Cov<-data.frame(cov1=c(4.336,5.615,4.770,4.655),cov2=c(1,0,0,1))
			 mle<-mle.vam(System&Time&Type~(GQR_ARAm(1.1,0.6|sqrt,3)|Weibull(0.1,2.1|0.6*cov1+(-0.7)*cov2))&(ARA1(0.9)+GQR_ARAm(1.3,0.8|log,2)),data=DataC,data.covariates = Cov)
			 theta<-c(0.11,2.5,1.15,0.62,0.92,1.25,0.75,0.65,-0.56)

			 lnL=-19.4945990891545
 dlnL=c(-206.714753771336,-6.82689872724557,-86.5113962999137,58.6590205517592,17.4686517107146,-28.1160439305364,12.3479412018079,-103.481756998376,-17.060552489827)
 d2lnL=matrix(c(-1818.18181818182,-46.205522372485,-1303.94531635136,919.348660858553,301.250316223058,-352.838567778391,182.230804898021,-1919.61597271251,-255.0959317257,-46.205522372485,-5.25170057388406,-56.4640109586639,46.4569611342398,15.2073077419991,-22.8395882079692,13.468952275987,-19.9452405296797,-8.09582203512599,-1303.94531635136,-56.4640109586639,-452.837947365564,323.283382515523,120.264062030744,-150.029357788677,74.4002758748625,-678.192600328426,-91.0803859824277,919.348660858553,46.4569611342398,323.283382515523,-285.512026132512,-86.3468653498148,83.400056974447,-57.1287622591366,478.335975277066,64.8090291013723,301.250316223058,15.2073077419991,120.264062030744,-86.3468653498148,-45.4752738368675,36.9912627317187,-20.0765210627799,154.711270435536,21.6881571086723,-352.838567778391,-22.8395882079692,-150.029357788677,83.400056974447,36.9912627317187,-34.5411711880569,25.5746498861266,-176.94480661365,-30.8959582022062,182.230804898021,13.468952275987,74.4002758748625,-57.1287622591366,-20.0765210627799,25.5746498861266,-14.512330314271,90.8052556393275,18.6791309672241,-1919.61597271251,-19.9452405296797,-678.192600328426,478.335975277066,154.711270435536,-176.94480661365,90.8052556393275,-1003.46078338031,-126.497958118298,-255.0959317257,-8.09582203512599,-91.0803859824277,64.8090291013723,21.6881571086723,-30.8959582022062,18.6791309672241,-126.497958118298,-28.060552489827),nrow=9,byrow=TRUE)
 C=-12.3714600692905
 dC=c(-4.24363873957824,-13.6103751907454,7.26004293245876,0.626338320171234,-8.38953273778471,2.15977488223833,3.84020105718946,-2.79864006881013)
 d2C=matrix(c(-4.08761269917479,-9.15929383105404,2.80523087133187,0.317682420416391,-7.27158052429994,2.8952684616044,1.98846788242069,-2.4134630610617,-9.15929383105404,-19.58313870284,8.43557460837913,5.92657696654463,-11.9923909476063,4.37714475266716,-0.595789462662424,-0.549291867754576,2.80523087133187,8.43557460837913,-39.0500515568676,-11.2971362801361,0.446501579833371,-7.74235024964267,0.506232798916784,0.678717820745405,0.317682420416391,5.92657696654463,-11.2971362801361,-23.6925636795104,3.37487513406574,-4.24302578662244,-0.83181929485292,0.444504989998187,-7.27158052429994,-11.9923909476063,0.446501579833371,3.37487513406574,-4.26553955031451,5.0447698657452,3.06925984071366,-3.22215810742276,2.8952684616044,4.37714475266716,-7.74235024964267,-4.24302578662244,5.0447698657452,-6.09776641393256,-1.8711642094494,3.00281437022592,1.98846788242069,-0.595789462662424,0.506232798916784,-0.83181929485292,3.06925984071366,-1.8711642094494,-3.36181186790202,2.92219115023161,-2.4134630610617,-0.549291867754576,0.678717820745405,0.444504989998187,-3.22215810742276,3.00281437022592,2.92219115023161,-5.14398244387502),nrow=8,byrow=TRUE)
 expect_that(logLik.mle.vam(mle,par0=theta) ,equals(lnL,tolerance=0.00000000000001))
 expect_that(contrast.mle.vam(mle,par0=theta) ,equals(C,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dlnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_gradient=TRUE) ,equals(dC,tolerance=0.00000000000001))
 expect_that(logLik.mle.vam(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2lnL,tolerance=0.00000000000001))
 expect_that(contrast(mle,par0=theta,with_value=FALSE,with_hessian=TRUE) ,equals(d2C,tolerance=0.00000000000001))
})
rcqls/VAM documentation built on Jan. 14, 2024, 9:07 p.m.