R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require("DoseFinding") Loading required package: DoseFinding > if(!require("multcomp")) + stop("need multcomp package to run this test") Loading required package: multcomp Loading required package: mvtnorm Loading required package: survival Loading required package: TH.data Loading required package: MASS Attaching package: 'TH.data' The following object is masked from 'package:MASS': geyser > > ######################################################################## > #### multContTest > # functions to sample random DF data > getDosSampSiz <- function(){ + # generate dose levels + mD <- runif(1, 0, 1500) + nD <- max(rpois(1, 5), 4) + p <- rgamma(nD, 3) + p <- cumsum(p/sum(p)) + doses <- signif(c(0, mD*p), 3) + # sample size allocations + totSS <- rpois(1, rexp(1, 1/250)) + totSS <- max(totSS, 50) + p <- rgamma(nD+1, 3);p <- p/sum(p) + n <- round(p*totSS) + n[n==0] <- rpois(sum(n==0), 1)+1 + list(doses=doses, n=n) + } > getDFdataSet <- function(doses, n){ + ll <- getDosSampSiz() + e0 <- rnorm(1, 0, 10) + eMax <- rgamma(1, abs(e0)*0.5, 0.5)*I(runif(1)<0.25) + if(eMax > 0){ sig <- eMax/runif(1, 0.5, 5)} + else { sig <- rgamma(1, abs(e0)*0.5, 0.5) } + dosVec <- rep(ll$doses, ll$n) + if(runif(1)<0.3){ + mnVec <- betaMod(dosVec, e0=e0, eMax=eMax, delta1=runif(1, 0.5, 5), + delta2=runif(1, 0.5, 5), scal=1.2*max(ll$doses)) + } else { + mnVec <- logistic(dosVec, e0 = e0, eMax = eMax, + ed50=runif(1, 0.05*max(ll$doses), 1.5*max(ll$doses)), + delta=runif(1, 0.5, max(ll$doses)/2)) + } + resp <- rnorm(sum(ll$n), mnVec, sig) + N <- sum(ll$n) + cov1 <- as.factor(rpois(N, 5)) + cov2 <- runif(N, 1, 100) + aa <- data.frame(x= dosVec, y=resp, cov1=cov1, cov2=cov2) + aa[sample(1:nrow(aa)),] + } > > #### simulate data and compare to output of glht of multcomp package and oldMCPMod function > set.seed(10) > dd <- getDFdataSet() > bet <- guesst(0.9*max(dd$x), p=0.8, "betaMod", scal = 1.2*max(dd$x), + dMax = 0.7*max(dd$x), Maxd = max(dd$x)) > sE <- guesst(c(0.5*max(dd$x), 0.7*max(dd$x)) , p=c(0.5, 0.9), "sigEmax") > models <- Mods(linear = NULL, betaMod = bet, sigEmax = sE, + doses = sort(unique(dd$x)), + addArgs=list(scal = 1.2*max(dd$x))) > obj <- MCTtest(x,y, dd, models=models, addCovars = ~cov1+cov2, pVal = T) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- lm(y~x+cov1+cov2, data=dd2) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = y ~ x + cov1 + cov2, data = dd2) Linear Hypotheses: Estimate Std. Error t value Pr(>t) linear <= 0 0.07260 0.05755 1.262 0.1756 betaMod <= 0 0.12179 0.05559 2.191 0.0287 * sigEmax <= 0 0.08677 0.05648 1.536 0.1114 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear betaMod sigEmax 0 -0.472 -0.492 -0.349 172 -0.391 -0.345 -0.489 285 -0.132 0.032 -0.240 516 0.262 0.799 0.503 761 0.733 0.007 0.574 Contrast Correlation: linear betaMod sigEmax linear 1.000 0.645 0.942 betaMod 0.645 1.000 0.741 sigEmax 0.942 0.741 1.000 Multiple Contrast Test: t-Stat adj-p betaMod 2.191 0.0287 sigEmax 1.536 0.1112 linear 1.262 0.1757 > > obj <- MCTtest(x,y, dd, models=models, addCovars = ~1, pVal = T) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- lm(y~x, data=dd2) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = y ~ x, data = dd2) Linear Hypotheses: Estimate Std. Error t value Pr(>t) linear <= 0 0.07985 0.05691 1.403 0.1404 betaMod <= 0 0.11839 0.05507 2.150 0.0314 * sigEmax <= 0 0.09407 0.05594 1.682 0.0853 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear betaMod sigEmax 0 -0.472 -0.487 -0.347 172 -0.389 -0.348 -0.490 285 -0.132 0.030 -0.239 516 0.255 0.800 0.500 761 0.738 0.005 0.576 Contrast Correlation: linear betaMod sigEmax linear 1.000 0.642 0.941 betaMod 0.642 1.000 0.740 sigEmax 0.941 0.740 1.000 Multiple Contrast Test: t-Stat adj-p betaMod 2.150 0.0317 sigEmax 1.682 0.0852 linear 1.403 0.1403 > > #### different model set > set.seed(10) > dd <- getDFdataSet() > mD <- max(dd$x) > lg1 <- guesst(c(0.3*mD, 0.4*mD), c(0.3, 0.9), "logistic") > lg2 <- guesst(c(0.3*mD, 0.4*mD), c(0.3, 0.5), "logistic") > expo <- guesst(c(0.9*mD), c(0.7), "exponential", Maxd=mD) > quad <- guesst(c(0.6*mD), c(1), "quadratic") > models <- Mods(linlog = NULL, logistic = rbind(lg1, lg2), + exponential = expo, quadratic = quad, + doses = sort(unique(dd$x)), addArgs=list(off = 0.2*max(dd$x))) > > obj <- MCTtest(x,y, dd, models=models, addCovars = ~cov1+cov2, pVal = T) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- lm(y~x+cov1+cov2, data=dd2) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = y ~ x + cov1 + cov2, data = dd2) Linear Hypotheses: Estimate Std. Error t value Pr(>t) linlog <= 0 0.09114 0.05971 1.526 0.1510 logistic1 <= 0 0.10293 0.05761 1.787 0.0940 . logistic2 <= 0 0.09583 0.05707 1.679 0.1154 exponential <= 0 0.02544 0.05514 0.461 0.5623 quadratic <= 0 0.11676 0.05980 1.952 0.0674 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linlog logistic1 logistic2 exponential quadratic 0 -0.624 -0.448 -0.457 -0.268 -0.665 172 -0.309 -0.590 -0.475 -0.347 -0.055 285 -0.049 0.101 -0.118 -0.191 0.199 516 0.366 0.500 0.499 -0.069 0.696 761 0.616 0.437 0.550 0.875 -0.174 Contrast Correlation: linlog logistic1 logistic2 exponential quadratic linlog 1.000 0.902 0.955 0.799 0.651 logistic1 0.902 1.000 0.949 0.671 0.660 logistic2 0.955 0.949 1.000 0.792 0.586 exponential 0.799 0.671 0.792 1.000 0.063 quadratic 0.651 0.660 0.586 0.063 1.000 Multiple Contrast Test: t-Stat adj-p quadratic 1.952 0.0675 logistic1 1.787 0.0941 logistic2 1.679 0.1154 linlog 1.526 0.1506 exponential 0.461 0.5625 > > obj <- MCTtest(x,y, dd, models=models, addCovars = ~1, pVal = T) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- lm(y~x, data=dd2) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = y ~ x, data = dd2) Linear Hypotheses: Estimate Std. Error t value Pr(>t) linlog <= 0 0.09484 0.05913 1.604 0.1320 logistic1 <= 0 0.10608 0.05706 1.859 0.0817 . logistic2 <= 0 0.10119 0.05654 1.790 0.0928 . exponential <= 0 0.03666 0.05440 0.674 0.4665 quadratic <= 0 0.10792 0.05938 1.817 0.0885 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linlog logistic1 logistic2 exponential quadratic 0 -0.623 -0.445 -0.455 -0.270 -0.664 172 -0.307 -0.592 -0.476 -0.342 -0.056 285 -0.050 0.101 -0.118 -0.189 0.197 516 0.358 0.497 0.496 -0.076 0.698 761 0.622 0.440 0.554 0.877 -0.175 Contrast Correlation: linlog logistic1 logistic2 exponential quadratic linlog 1.000 0.901 0.955 0.800 0.648 logistic1 0.901 1.000 0.949 0.670 0.658 logistic2 0.955 0.949 1.000 0.792 0.584 exponential 0.800 0.670 0.792 1.000 0.062 quadratic 0.648 0.658 0.584 0.062 1.000 Multiple Contrast Test: t-Stat adj-p logistic1 1.859 0.0816 quadratic 1.817 0.0887 logistic2 1.790 0.0933 linlog 1.604 0.1316 exponential 0.674 0.4667 > > #### contrast matrix handed over > set.seed(23) > dd <- getDFdataSet() > mD <- max(dd$x) > lg1 <- guesst(c(0.3*mD, 0.4*mD), c(0.3, 0.9), "logistic") > lg2 <- guesst(c(0.3*mD, 0.4*mD), c(0.3, 0.5), "logistic") > expo <- guesst(c(0.9*mD), c(0.7), "exponential", Maxd=mD) > quad <- guesst(c(0.6*mD), c(1), "quadratic") > models <- Mods(linlog = NULL, logistic = rbind(lg1, lg2), + exponential = expo, quadratic = quad, + doses = dd$x, addArgs=list(off = 0.2*max(dd$x))) > > obj <- MCTtest(x,y, dd, models=models, addCovars = ~cov1+cov2, pVal = T) > contMat <- obj$contMat > obj <- MCTtest(x,y, dd, models=models, addCovars = ~1, pVal = T, contMat = contMat) > dd2 <- dd > dd2$x <- as.factor(dd2$x) > fit <- lm(y~x, data=dd2) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: lm(formula = y ~ x, data = dd2) Linear Hypotheses: Estimate Std. Error t value Pr(>t) linlog <= 0 -0.2995 1.0212 -0.293 0.821 logistic1 <= 0 -0.1795 0.9837 -0.182 0.785 logistic2 <= 0 -0.1401 0.9877 -0.142 0.771 exponential <= 0 -0.2892 0.9568 -0.302 0.823 quadratic <= 0 -0.1973 1.2830 -0.154 0.776 (Adjusted p values reported -- single-step method) > obj Multiple Contrast Test Contrasts: linlog logistic1 logistic2 exponential quadratic 0 -0.464 -0.268 -0.289 -0.207 -0.751 124 -0.436 -0.627 -0.556 -0.450 -0.241 272 0.013 -0.006 -0.083 -0.187 0.476 560 0.127 0.196 0.172 -0.004 0.353 865 0.761 0.705 0.755 0.848 0.164 Contrast Correlation: linlog logistic1 logistic2 exponential quadratic linlog 1.000 0.950 0.963 0.890 0.665 logistic1 0.950 1.000 0.992 0.879 0.591 logistic2 0.963 0.992 1.000 0.922 0.542 exponential 0.890 0.879 0.922 1.000 0.252 quadratic 0.665 0.591 0.542 0.252 1.000 Multiple Contrast Test: t-Stat adj-p logistic2 -0.142 0.771 quadratic -0.154 0.776 logistic1 -0.182 0.785 linlog -0.293 0.821 exponential -0.302 0.823 > > ######################################################################## > #### some binary test cases > getDFdataSet.bin <- function(doses, n){ + ll <- getDosSampSiz() + ll$n <- ll$n+10 + e0 <- rnorm(1, 0, sqrt(3.28)) + eMax <- rnorm(1, 0, 5) + dosVec <- rep(ll$doses, ll$n) + if(runif(1)<0.3){ + mn <- betaMod(dosVec, e0 = e0, eMax = eMax, delta1=runif(1, 0.5, 5), + delta2=runif(1, 0.5, 5), scal=1.2*max(ll$doses)) + } else { + mn <- logistic(dosVec, e0 = e0, + eMax = eMax, ed50=runif(1, 0.05*max(ll$doses), 1.5*max(ll$doses)), + delta=runif(1, 0.5, max(ll$doses)/2)) + } + resp <- rbinom(length(ll$n), ll$n, 1/(1+exp(-mn))) + aa <- data.frame(dose = ll$doses, resp = resp) + aa <- data.frame(x= aa$dose, y=aa$resp/ll$n, n=ll$n) + aa[sample(1:nrow(aa)),] + } > > set.seed(1909) > dd <- getDFdataSet.bin() > bet <- guesst(0.9*max(dd$x), p=0.8, "betaMod", scal = 1.2*max(dd$x), dMax = 0.7*max(dd$x), + Maxd = max(dd$x)) > sE <- guesst(c(0.5*max(dd$x), 0.7*max(dd$x)) , p=c(0.5, 0.9), "sigEmax") > models <- Mods(linear = NULL, betaMod = bet, sigEmax = sE, + doses = sort(unique(dd$x)), addArgs=list(scal = 1.2*max(dd$x))) > logReg <- glm(y~as.factor(x)-1, family=binomial, data=dd, weights = n) > dePar <- coef(logReg) > vCov <- vcov(logReg) > dose <- sort(unique(dd$x)) > obj <- MCTtest(dose, dePar, S=vCov, models=models, type="general", + df=Inf, pVal = T) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- glm(y~x-1, family = binomial, data=dd2, weights = n) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: glm(formula = y ~ x - 1, family = binomial, data = dd2, weights = n) Linear Hypotheses: Estimate Std. Error z value Pr(>z) linear <= 0 0.19645 0.33561 0.585 0.395 betaMod <= 0 0.03813 0.31565 0.121 0.587 sigEmax <= 0 0.30472 0.34592 0.881 0.283 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear betaMod sigEmax 0 -0.794 -0.793 -0.774 190 -0.064 -0.022 -0.192 344 0.197 0.577 0.317 404 0.098 0.190 0.163 616 0.563 0.048 0.487 Contrast Correlation: linear betaMod sigEmax linear 1.000 0.723 0.962 betaMod 0.723 1.000 0.793 sigEmax 0.962 0.793 1.000 Multiple Contrast Test: t-Stat adj-p sigEmax 0.881 0.283 linear 0.585 0.395 betaMod 0.121 0.587 > > set.seed(1997) > dd <- getDFdataSet.bin() > bet <- guesst(0.9*max(dd$x), p=0.8, "betaMod", scal = 1.2*max(dd$x), + dMax = 0.7*max(dd$x), Maxd = max(dd$x)) > sE <- guesst(c(0.5*max(dd$x), 0.7*max(dd$x)) , p=c(0.5, 0.9), "sigEmax") > models <- Mods(linear = NULL, betaMod = bet, sigEmax = sE,direction = "decreasing", + addArgs=list(scal = 1.2*max(dd$x)), doses = sort(unique(dd$x))) > logReg <- glm(y~as.factor(x)-1, family=binomial, data=dd, weights = n) > dePar <- coef(logReg) > vCov <- vcov(logReg) > dose <- sort(unique(dd$x)) > obj <- MCTtest(dose, dePar, S=vCov, models=models, type = "general", + pVal = TRUE, df=Inf) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- glm(y~x-1, family = binomial, data=dd2, weights = n) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: glm(formula = y ~ x - 1, family = binomial, data = dd2, weights = n) Linear Hypotheses: Estimate Std. Error z value Pr(>z) linear <= 0 0.2547 0.1877 1.357 0.1500 betaMod <= 0 0.2367 0.2371 0.999 0.2555 sigEmax <= 0 0.3698 0.1889 1.958 0.0476 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear betaMod sigEmax 0 0.418 0.634 0.310 68.5 0.114 0.150 0.139 85.9 0.387 0.362 0.551 166 0.027 -0.572 -0.033 197 -0.015 -0.158 -0.054 278 -0.129 -0.264 -0.172 332 -0.803 -0.152 -0.740 Contrast Correlation: linear betaMod sigEmax linear 1.000 0.635 0.967 betaMod 0.635 1.000 0.677 sigEmax 0.967 0.677 1.000 Multiple Contrast Test: t-Stat adj-p sigEmax 1.958 0.0477 linear 1.357 0.1501 betaMod 0.999 0.2556 > > set.seed(1) > dd <- getDFdataSet.bin() > bet <- guesst(0.9*max(dd$x), p=0.8, "betaMod", scal = 1.2*max(dd$x), + dMax = 0.7*max(dd$x), Maxd = max(dd$x)) > sE <- guesst(c(0.5*max(dd$x), 0.7*max(dd$x)) , p=c(0.5, 0.9), "sigEmax") > models <- Mods(linear = NULL, betaMod = bet, sigEmax = sE, + doses = sort(unique(dd$x)), addArgs=list(scal = 1.2*max(dd$x))) > logReg <- glm(y~as.factor(x)-1, family=binomial, data=dd, weights = n) > dePar <- coef(logReg) > vCov <- vcov(logReg) > dose <- sort(unique(dd$x)) > obj <- MCTtest(dose, dePar, S=vCov, models=models, type = "general", + pVal = T, df=Inf) > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- glm(y~x-1, family = binomial, data=dd2, weights = n) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: glm(formula = y ~ x - 1, family = binomial, data = dd2, weights = n) Linear Hypotheses: Estimate Std. Error z value Pr(>z) linear <= 0 0.3391 0.3218 1.054 0.226 betaMod <= 0 0.3556 0.3124 1.138 0.201 sigEmax <= 0 0.2935 0.3497 0.839 0.299 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear betaMod sigEmax 0 -0.721 -0.796 -0.547 79.1 -0.129 -0.102 -0.312 117 0.005 0.106 -0.137 208 0.191 0.519 0.288 398 0.654 0.274 0.708 Contrast Correlation: linear betaMod sigEmax linear 1.000 0.785 0.954 betaMod 0.785 1.000 0.742 sigEmax 0.954 0.742 1.000 Multiple Contrast Test: t-Stat adj-p betaMod 1.138 0.201 linear 1.054 0.226 sigEmax 0.839 0.299 > > ## one-dimensional test > set.seed(1) > dd <- getDFdataSet.bin() > model <- Mods(linear = NULL, doses=sort(unique(dd$x)), addArgs=list(scal = 1.2*max(dd$x))) > logReg <- glm(y~as.factor(x)-1, family=binomial, data=dd, weights = n) > dePar <- coef(logReg) > vCov <- vcov(logReg) > dose <- sort(unique(dd$x)) > obj <- MCTtest(dose, dePar, S=vCov, models=model, type = "general", + pVal = T, df=Inf) Warning message: In MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol) : Warning from mvtnorm::pmvt: univariate: using pnorm. > dd2 <- dd;dd2$x <- as.factor(dd$x) > fit <- glm(y~x-1, family = binomial, data=dd2, weights = n) > mcp <- glht(fit, linfct = mcp(x = t(obj$contMat)), alternative = "greater") > summary(mcp) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: glm(formula = y ~ x - 1, family = binomial, data = dd2, weights = n) Linear Hypotheses: Estimate Std. Error z value Pr(>z) linear <= 0 0.3391 0.3218 1.054 0.146 (Adjusted p values reported -- single-step method) > print(obj, digits = 3) Multiple Contrast Test Contrasts: linear 0 -0.721 79.1 -0.129 117 0.005 208 0.191 398 0.654 Contrast Correlation: linear linear 1 Multiple Contrast Test: t-Stat adj-p linear 1.054 0.146 > > > ######################################################################## > ## unordered values in MCTtest > ## placebo-adjusted scale > ## two blocks below should give equal results > data(IBScovars) > modlist <- Mods(emax = 0.05, linear = NULL, logistic = c(0.5, 0.1), + linInt = c(0, 1, 1, 1), doses = c(0, 1, 2, 3, 4)) > ancMod <- lm(resp~factor(dose)+gender, data=IBScovars) > drEst <- coef(ancMod)[2:5] > vc <- vcov(ancMod)[2:5, 2:5] > doses <- 1:4 > fitMod(doses, drEst, S=vc, model = "sigEmax", placAdj=TRUE, type = "general") Message: Need bounds in "bnds" for nonlinear models, using default bounds from "defBnds". Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model eMax ed50 h 0.4714 0.4880 0.5000 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.07825 > MCTtest(doses, drEst, S = vc, models = modlist, placAdj = TRUE, + type = "general", df = Inf) Multiple Contrast Test Contrasts: emax linear logistic linInt 1 0.462 -0.429 0.511 -0.672 2 0.505 0.002 0.507 0.437 3 0.504 0.400 0.485 0.418 4 0.527 0.810 0.497 0.428 Contrast Correlation: emax linear logistic linInt emax 1.000 0.715 1.000 0.617 linear 0.715 1.000 0.698 0.863 logistic 1.000 0.698 1.000 0.598 linInt 0.617 0.863 0.598 1.000 Multiple Contrast Test: t-Stat adj-p emax 3.178 0.00193 logistic 3.164 0.00194 linear 2.640 0.00968 linInt 2.247 0.02721 > > > ord <- c(3,4,1,2) > drEst2 <- drEst[ord] > vc2 <- vc[ord,ord] > doses2 <- doses[ord] > fitMod(doses2, drEst2, S=vc2, model = "sigEmax", placAdj=TRUE, type = "general") Message: Need bounds in "bnds" for nonlinear models, using default bounds from "defBnds". Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model eMax ed50 h 0.4714 0.4880 0.5000 Fitted to: 3 4 1 2 0.3502 0.3480 0.2846 0.2965 Generalized residual sum of squares: 0.07825 > MCTtest(doses2, drEst2, S = vc2, models = modlist, placAdj = TRUE, + type = "general", df = Inf) Multiple Contrast Test Contrasts: emax linear logistic linInt 1 0.462 -0.429 0.511 -0.672 2 0.505 0.002 0.507 0.437 3 0.504 0.400 0.485 0.418 4 0.527 0.810 0.497 0.428 Contrast Correlation: emax linear logistic linInt emax 1.000 0.715 1.000 0.617 linear 0.715 1.000 0.698 0.863 logistic 1.000 0.698 1.000 0.598 linInt 0.617 0.863 0.598 1.000 Multiple Contrast Test: t-Stat adj-p emax 3.178 0.00213 logistic 3.164 0.00174 linear 2.640 0.00941 linInt 2.247 0.02663 > > ## unadjusted scale > ## two blocks below should give equal results > ancMod <- lm(resp~factor(dose)-1, data=IBScovars) > drEst <- coef(ancMod) > vc <- vcov(ancMod) > doses <- 0:4 > fitMod(doses, drEst, S=vc, model = "sigEmax", type = "general") Message: Need bounds in "bnds" for nonlinear models, using default bounds from "defBnds". Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model e0 eMax ed50 h 0.2171 0.4714 0.4860 0.5000 Fitted to: 0 1 2 3 4 0.2169 0.5016 0.5138 0.5677 0.5648 Generalized residual sum of squares: 0.07892 > MCTtest(doses, drEst, S = vc, models = modlist, type = "general", df = Inf) Multiple Contrast Test Contrasts: emax linear logistic linInt 0 -0.894 -0.617 -0.894 -0.522 1 0.207 -0.338 0.228 -0.573 2 0.226 0.002 0.227 0.373 3 0.227 0.315 0.218 0.358 4 0.235 0.637 0.221 0.363 Contrast Correlation: emax linear logistic linInt emax 1.000 0.715 1.000 0.618 linear 0.715 1.000 0.698 0.863 logistic 1.000 0.698 1.000 0.599 linInt 0.618 0.863 0.599 1.000 Multiple Contrast Test: t-Stat adj-p emax 3.186 0.00184 logistic 3.172 0.00177 linear 2.645 0.00956 linInt 2.254 0.02649 > > > ord <- c(3,4,1,2,5) > drEst2 <- drEst[ord] > vc2 <- vc[ord,ord] > doses2 <- doses[ord] > fitMod(doses2, drEst2, S=vc2, model = "sigEmax", type = "general") Message: Need bounds in "bnds" for nonlinear models, using default bounds from "defBnds". Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model e0 eMax ed50 h 0.2171 0.4714 0.4860 0.5000 Fitted to: 2 3 0 1 4 0.5138 0.5677 0.2169 0.5016 0.5648 Generalized residual sum of squares: 0.07892 > MCTtest(doses2, drEst2, S = vc2, models = modlist, type = "general", df = Inf) Multiple Contrast Test Contrasts: emax linear logistic linInt 0 -0.894 -0.617 -0.894 -0.522 1 0.207 -0.338 0.228 -0.573 2 0.226 0.002 0.227 0.373 3 0.227 0.315 0.218 0.358 4 0.235 0.637 0.221 0.363 Contrast Correlation: emax linear logistic linInt emax 1.000 0.715 1.000 0.618 linear 0.715 1.000 0.698 0.863 logistic 1.000 0.698 1.000 0.599 linInt 0.618 0.863 0.599 1.000 Multiple Contrast Test: t-Stat adj-p emax 3.186 0.00163 logistic 3.172 0.00192 linear 2.645 0.00981 linInt 2.254 0.02699 > > ######################################################################## > ## catch cases where mvtnorm does not calculate result due to non-psd > ## covariance matrix > doses<-c(0,10,20,40) > exm1<-0.15 > exm2<-c(0.05,5) > expo<-0.2 > quad<--0.6 > beta<-c(0.05,4) > data.sim <- structure(list(X = structure(1:4, .Label = c("0", "10", "20", "40"), class = "factor"), + dose = c(0L, 10L, 20L, 40L), + Estimate = c(0.266942236, 3.792703657, 14.69084734, 17.71179102), + Cov1 = c(3.685607913, 0.595285049, 0.651289991, 0.742901538), + Cov2 = c(0.595285049, 3.31255546, 0.47843908, 0.545737127), + Cov3 = c(0.651289991, 0.47843908, 3.398708786, 0.597080557), + Cov4 = c(0.742901538, 0.545737127, 0.597080557, 3.556324648)), + class = "data.frame", row.names = c(NA, -4L)) > mu<-data.sim[,3] > S<-data.matrix(data.sim[,4:7],rownames.force = NA) > models2<-Mods(doses=doses, emax=exm1,sigEmax=exm2,linear=NULL,exponential=expo,quadratic=quad,betaMod=beta) > tst <- MCTtest(dose=doses,resp=mu,models = models2,S=S,type="general") Warning messages: 1: In MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol) : Warning from mvtnorm::pmvt: Covariance matrix not positive semidefinite. 2: In MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol) : Setting calculated p-value to NA > ## p-value of linear model should be NA > is.na(attr(tst$tStat, "pVal")[3]) [1] TRUE > > proc.time() user system elapsed 3.64 0.34 3.95