R Under development (unstable) (2023-10-30 r85440 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.

> ## commented out for time (and dependency reasons)
> require("DoseFinding")
Loading required package: DoseFinding
Loading required package: ggplot2
Loading required package: lattice
Loading required package: mvtnorm
> if(!(require("quadprog") & require("Rsolnp")))
+   stop("need packages quadprog and Rsolnp to run these tests")
Loading required package: quadprog
Loading required package: Rsolnp
> 
> ## ## calculation of optimal contrast by enumerating all active sets
> ## allActiveSets <- function(S, mu, mult){
> ##   k <- length(mu)
> ##   CC <- cbind(-1, diag(k - 1))
> ##   SPa <- CC %*% S %*% t(CC)
> ##   muPa <- as.numeric(CC %*% mu)
> ##   ## generate all possible active sets
> ##   mat <- matrix(nrow = 2^(k-1), ncol = (k-1))
> ##   for(i in 1:(k-1))
> ##     mat[,i] <- rep(rep(c(FALSE,TRUE), each=2^(i-1)), 2^((k-1)-i))
> ##   val <- numeric(2^(k-1))
> ##   feasible <- logical(2^(k-1))
> ##   cont <- matrix(nrow = 2^(k-1), ncol = (k-1))
> ##   for(i in 1:(2^(k-1))){
> ##     nonzero <- mat[i,]
> ##     if(sum(nonzero) > 0){
> ##       cont[i,!nonzero] <- 0
> ##       cont[i,nonzero] <- solve(SPa[nonzero, nonzero]) %*% muPa[nonzero]
> ##       feasible[i] <- all(mult*cont[i,] >= 0)
> ##       contrast <- c(-sum(cont[i,]), cont[i,])
> ##       val[i] <- as.numeric(t(contrast)%*%mu/sqrt(t(contrast)%*%S%*%contrast))
> ##     }
> ##   }
> ##   if(!any(feasible))
> ##     return(rep(NA, k))
> ##   mm <- max(val[which(feasible)])
> ##   c(-sum(cont[val == mm,]), cont[val == mm,])  
> ## }
> 
> 
> ## ## helper functions
> ## getStand <- function(x)
> ##   x/sqrt(sum(x^2))
> ## getNCP <- function(cont, mu, S)
> ##   as.numeric(t(cont)%*%mu/sqrt(t(cont)%*%S%*%cont))
> 
> ## set.seed(1)
> ## ncp1 <- ncp2 <- ncp3 <- ncp4 <- ncp5 <- numeric(1000)
> ## for(i in 1:1000){
> ##   ## simulate mean and covariance matrix
> ##   kk <- round(runif(1, 4, 10))
> ##   A <- matrix(runif(kk^2,-1,1),kk,kk)
> ##   S <- crossprod(A)+diag(kk)
> ##   mult <- sign(rnorm(1))
> ##   mu <- mult*sort(rnorm(kk, 1:kk, 1))
>   
> ##   ## unconstrained solution
> ##   ones <- rep(1,kk)
> ##   unConst <- solve(S)%*%(mu - c(t(mu)%*%solve(S)%*%ones/(t(ones)%*%solve(S)%*%ones)))
> ##   cont1 <- getStand(unConst)
>   
> ##   ## function from DoseFinding package
> ##   cont2 <- DoseFinding:::constOptC(mu, solve(S), placAdj=FALSE,
> ##                                    ifelse(mult == 1, "increasing", "decreasing"))
>   
> ##   ## alternative solution using quadratic programming
> ##   D <- S
> ##   d <- rep(0,kk)
> ##   tA <- rbind(rep(1, kk), 
> ##               mu, 
> ##               mult*diag(kk)*c(-1,rep(1,kk-1)))
> ##   A <- t(tA)
> ##   bvec <- c(0,1,rep(0,kk))
> ##   rr <- solve.QP(D, d, A, bvec, meq=2)
> ##   cont3 <- getStand(rr$solution)
>   
> ##   ## using solnp
> ##   LB <- rep(0, kk-1)
> ##   UB <- rep(20, kk-1)
> ##   strt <- rep(1, kk-1)
> ##   mgetNCP <- function(x, ...){
> ##     cont <- c(-sum(x), x)
> ##     -getNCP(cont, ...)
> ##   }
> ##   res <- solnp(strt, mgetNCP, mu=mu, S=S,
> ##                LB=LB, UB=UB,
> ##                control = list(trace = 0))
> ##   out <- c(-sum(res$pars), res$pars)
> ##   cont4 <- getStand(out)
> 
> ##   ## using
> ##   cont5 <- allActiveSets(S=S, mu=mu, mult=mult)
>     
> ##   ## compare optimized non-centrality parameters
> ##   ncp1[i] <- getNCP(cont1, mu, S)
> ##   ncp2[i] <- getNCP(cont2, mu, S)
> ##   ncp3[i] <- getNCP(cont3, mu, S)
> ##   ncp4[i] <- getNCP(cont4, mu, S)
> ##   ncp5[i] <- getNCP(cont5, mu, S)  
> ## }
> ## sapply(list(ncp1, ncp2, ncp3, ncp4, ncp5), quantile)
> 
> ## ## tests whether constant shapes (possible with linInt) are handled correctly
> ## data(biom)
> ## ## define shapes for which to calculate optimal contrasts
> ## modlist <- Mods(emax = 0.05, linear = NULL, logistic = c(0.5, 0.1),
> ##                 linInt = rbind(c(0, 0, 0, 1), c(0, 1, 1, 1)),
> ##                 doses = c(0, 0.05, 0.2, 0.6, 1), placEff = 1)
> 
> ## optContr(modlist, w=1, doses=c(0.05), placAdj=TRUE, type = "u")
> ## optContr(modlist, w=1, doses=c(0.05), placAdj=TRUE, type = "c")
> ## optContr(modlist, w=1, doses=c(0.05,0.5), placAdj=TRUE, type = "u")
> ## optContr(modlist, w=1, doses=c(0.05,0.5), placAdj=TRUE, type = "c")
> ## optContr(modlist, w=1, doses=c(0,0.05), placAdj=FALSE, type = "u")
> ## optContr(modlist, w=1, doses=c(0,0.05), placAdj=FALSE, type = "c")
> ## optContr(modlist, w=1, doses=c(0,0.05,0.5), placAdj=FALSE, type = "u")
> ## optContr(modlist, w=1, doses=c(0,0.05,0.5), placAdj=FALSE, type = "c")
> 
> 
> ## modlist2 <- Mods(linInt = rbind(c(0, 1, 1, 1), c(0,0,0,1)),
> ##                  doses = c(0, 0.05, 0.2, 0.6, 1), placEff = 1)
> 
> ## ## all of these should throw an error
> ## optContr(modlist2, w=1, doses=c(0.05), placAdj=TRUE, type = "u")
> ## optContr(modlist2, w=1, doses=c(0.05), placAdj=TRUE, type = "c")
> ## optContr(modlist2, w=1, doses=c(0,0.05), placAdj=FALSE, type = "u")
> ## optContr(modlist2, w=1, doses=c(0,0.05), placAdj=FALSE, type = "c")
> ## ## these should work
> ## optContr(modlist2, w=1, doses=c(0.05,0.5), placAdj=TRUE, type = "u")
> ## optContr(modlist2, w=1, doses=c(0.05,0.5), placAdj=TRUE, type = "c")
> ## optContr(modlist2, w=1, doses=c(0,0.05,0.5), placAdj=FALSE, type = "u")
> ## optContr(modlist2, w=1, doses=c(0,0.05,0.5), placAdj=FALSE, type = "c")
> 
> proc.time()
   user  system elapsed 
   0.68    0.07    0.75