tests/tinytest/test-migpdCoefs.R

In texmex: Statistical Modelling of Extreme Values

library("MASS")
liver <- liver
liver$ndose <- as.numeric(liver$dose)

ralt <- resid(rlm(log(ALT.M) ~ log(ALT.B) + ndose, data=liver))
rast <- resid(rlm(log(AST.M) ~ log(AST.B) + ndose, data=liver))
ralp <- resid(rlm(log(ALP.M) ~ log(ALP.B) + ndose, data=liver))
rtbl <- resid(rlm(log(TBL.M) ~ log(TBL.B) + ndose, data=liver))

rliver <- data.frame(alt=ralt, ast=rast, alp=ralp, tbl=rtbl, ndose=liver$ndose)

Dmod <- migpd(rliver[rliver$ndose == 4, 1:4], mqu=.7) # Model for dose D

oldALTco <- coef(Dmod)[3:4, 1]

altgpd <- evm(alt, qu=.7, xi = ~ ndose, data=rliver)
astgpd <- evm(ast, qu=.7, xi = ~ ndose, data=rliver)

altco <- c(coef(altgpd)[1], coef(altgpd)[2] + 4 * coef(altgpd)[3])
astco <- c(coef(astgpd)[1], coef(astgpd)[2] + 4 * coef(astgpd)[3])

# Change one set of coefficients
lmod <- migpdCoefs(Dmod, which="alt", list(altco))

newALTco <- coef(lmod)[3:4, 1]
newALTco[1] <- log(newALTco[1]) # For comparison with altco
oldALTco[1] <- log(oldALTco[1])

expect_equivalent(altco, newALTco,
                  info="migpdCoefs: change one set of coefficients")
expect_true(all(newALTco!=oldALTco),
            info="migpdCoefs: change one set of coefficients")

# Change 2 sets of coefficients at once

lmod <- migpdCoefs(Dmod, which=c("alt", "ast"), coefs=list(altco, astco))

newCo <- coef(lmod)[3:4, 1:2]
oldCo <- coef(Dmod)[3:4, 1:2]

newCo[1,] <- log(newCo[1,])
oldCo[1,] <- log(oldCo[1,])

expect_equal(c(cbind(altco, astco)), c(newCo),
             info="migpdCoefs: change two sets of coefficients at once")
expect_true(all(newCo!=oldCo),
            info="migpdCoefs: change two set of  coefficients at once")