Nothing
R/LRcontrast.R
In LRcontrast: Dose Response Signal Detection under Model Uncertainty
#SIMULATING QUANTILES:
#define variables:
#not changeable for the user:
env1<-new.env() #environment in which scal, off, i is saved
SCALE1<-c(1,1,1,1) #scaling of the parameters in fct optim (parscale)
LOW1<-c(0,0,0,0) #lower bounds for (phi1, phi2, theta1, theta2)
UP1<-c(pi,pi/2,6,6) #upper bounds for (phi1, phi2, theta) (LOW and UP needed in nlminb)
INTERVAL<-c(0, pi) #interval for phi (needed in optimize)
#starting values for nlminb in qLRcontrast and pLRcontrast:
#for the models linear, linlog, quadratic, emax, exponential:
initVal1<-matrix(nrow=3,ncol=3) #columns of the matrix are different starting values for (phi1, phi2, theta) in nlminb
initVal1[3,]<-rep(2,3)
initVal1[2,]<-rep(1,3)
initVal1[1,]<-c(0,pi/2,pi)
#for the models logistic and betaMod
initVal2<-matrix(nrow=4,ncol=8)
initVal2[4,]<-rep(2,8)
initVal2[3,]<-rep(2,8)
initVal2[2,]<-rep(1,8)
initVal2[1,]<-c(0,0.5,1,1.5,2,2.5,3,pi)
#for the model sigEmax
initVal3<-matrix(nrow=4,ncol=4)
initVal3[4,]<-rep(2,4)
initVal3[3,]<-rep(2,4)
initVal3[2,]<-rep(1,4)
initVal3[1,]<-c(0,1,2,pi)
#define the stochastic process:
#factor (-1), because nlminb minimizes functions instead of maximizing ( max{f(x)} = -min{-f(x)} )):
#parametrization with polar coordinates:
#here: beta0 = cos(phi), sigma = sin(phi)
WS1<-function(phi){return(sum(sqrt(env1$weight)*(cos(phi)*env1$Z[,env1$i]+sqrt(2)*sin(phi)*env1$Y[,env1$i]))/sqrt(sum(env1$weight*(cos(phi)^2+2*sin(phi)^2))))}
#here: beta0 = cos(phi1)*cos(phi2), beta1 = sin(phi2), sigma = sin(phi1)*cos(phi2)
WS2_linear<-function(x){return(.C("WS2linear", phi1=as.double(x[1]), phi2=as.double(x[2]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_emax<-function(x){return(.C("WS2emax", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_exponential<-function(x){return(.C("WS2exponential", phi1=as.double(x[1]), phi2=as.double(x[2]), delta=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_linlog<-function(x){return(.C("WS2linlog", phi1=as.double(x[1]), phi2=as.double(x[2]), off=as.double(env1$off), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_sigEmax<-function(x){return(.C("WS2sigEmax", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), h=as.double(x[4]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_quadratic<-function(x){return(.C("WS2quadratic", phi1=as.double(x[1]), phi2=as.double(x[2]), b=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_betaMod<-function(x){return(.C("WS2betaMod", phi1=as.double(x[1]), phi2=as.double(x[2]), delta1=as.double(x[3]), delta2=as.double(x[4]), scal=as.double(env1$scal), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_logistic<-function(x){return(.C("WS2logistic", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), delta=as.double(x[4]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
#maximizing the functions:
#maximize WS1 over INTERVAL=[0,2*pi] (M1) (equal for all models)
M1max<-function(x){return(max(optimize(WS1, interval=INTERVAL, maximum=TRUE)$objective, 0)^2)}
#maximize WS2 for the different models over M2. Maximize with nlminb with different starting values and take the maximum.
MAX_linear<-function(x){max(max(0,(-1)*nlminb(initVal1[1:2,1], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,2], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,3], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_emax<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_exponential<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_linlog<-function(x){max(max(0,(-1)*nlminb(initVal1[1:2,1], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,2], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,3], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_quadratic<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_betaMod<-function(x){max(max(0,(-1)*nlminb(initVal2[,1], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,2], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,3], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,4], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,5], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,6], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,7], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,8], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_logistic<-function(x){max(max(0,(-1)*nlminb(initVal2[,1], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,2], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,3], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,4], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,5], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,6], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,7], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,8], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_sigEmax<-function(x){max(max(0,(-1)*nlminb(initVal3[,1], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,2], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,3], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,4], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
#with that function the user can simulate the quantiles:
qLRcontrast<-function(dose, probs, models, weight = rep(1/length(dose), length(dose)), off = 0.01 * max(dose), scal = 1.2 * max(dose), nsim = 10000, info = TRUE)
{
env1$weight<-weight;
if(sum(env1$weight)!=1){stop("The sum of 'weight' needs to be 1")};
if(min(env1$weight)<=0){stop("Only positive values in 'weight' allowed")};
env1$dose<-dose;
if(min(env1$dose) < 0){stop("Only dose-levels >= 0 allowed")};
if(length(env1$dose)!=length(env1$weight)){stop("Lengths of the vectors 'dose' and 'weight' need to be equal")};
env1$probs<-sort(probs);
if(env1$probs[1]<0 | 1<max(env1$probs)){stop("'probs' needs to be a numeric vector of probabilities with values in [0,1]")};
env1$off<-off;
if(env1$off <= 0){stop("'off' parameter needs to be positive")};
env1$scal<-scal;
if(env1$scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
if(nsim<=0 | nsim%%1!=0){stop("'nsim' needs to be a positive integer")};
env1$Z<-matrix(data=rnorm(length(dose)*nsim),nrow=length(dose));
env1$Y<-matrix(data=rnorm(length(dose)*nsim),nrow=length(dose));
env1$L<-matrix(nrow=length(models),ncol=nsim);
L_max<-c();
QUANTILE<-matrix(nrow=length(models)+1, ncol=length(env1$probs)); colnames(QUANTILE)<-env1$probs; rownames(QUANTILE)<-c(models, "max");
maxFkt<-list()
for(j in 1:length(models))
{
switch(models[j],
linear={maxFkt[[j]]<-MAX_linear},
emax={maxFkt[[j]]<-MAX_emax},
exponential={maxFkt[[j]]<-MAX_exponential},
linlog={maxFkt[[j]]<-MAX_linlog},
sigEmax={maxFkt[[j]]<-MAX_sigEmax},
quadratic={maxFkt[[j]]<-MAX_quadratic},
betaMod={maxFkt[[j]]<-MAX_betaMod},
logistic={maxFkt[[j]]<-MAX_logistic},
stop("The vector 'models' contains an incorrect entry")
)
};
start<-Sys.time();
if(info==TRUE){cat(paste("% done \t", "nsim done\t", "secs passed\t", "secs remaining ~", collapse=""), sep="", fill=TRUE)};
for(j in 1:nsim){env1$i<-j; M1_max<-M1max(1);
for(k in 1:length(models)){env1$L[k,j]<-maxFkt[[k]](1)-M1_max};
L_max[j]<-max(env1$L[,j])
if(info==TRUE){
passed <- round(difftime(Sys.time(),start,units="secs")[[1]],0);
percent<- round(j/nsim*100,2);
if(j%%10==0){env1$remaining<- round(passed/j*(nsim-j),0)};
cat(paste(rep.int("\b", 60),collapse=""),sep="");
cat(paste(percent, "\t\t", j, "\t\t", passed, "\t\t", env1$remaining, collapse=""), sep="");
flush.console()}
};
for(j in 1:length(models)){QUANTILE[j,]<-quantile(env1$L[j,],env1$probs)};
QUANTILE[length(models)+1,]<-quantile(L_max,env1$probs);
cat("\n");
return(QUANTILE)}
#test statistic (2.9):
sLRcontrast<-function(dose, resp, models, off = 0.01 * max(dose), scal= 1.2 * max(dose))
{
if(length(dose)!=length(resp)){stop("The length of the vectors 'dose' and 'resp' needs to be equal")};
if(min(dose) < 0){stop("Only dose-levels >= 0 allowed")};
if(off <= 0){stop("'off' parameter needs to be positive")};
if(scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
result <- matrix(nrow=length(models)+1, ncol= 1); rownames(result)<-c(models, "max"); colnames(result)<-c("statistic");
RSS <- c();
addArgs <- list(); addArgs$off <- off; addArgs$scal <- scal;
bnds <- DoseFinding::defBnds(mD = max(dose));
minT0<-sum((resp-1/length(resp)*sum(resp))^2);
for(j in 1:length(models))
{
switch(models[j],
linear = {fit <- DoseFinding::fitMod(dose, resp, model = "linear", addArgs = addArgs, bnds = bnds$linear);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
emax = {fit <- DoseFinding::fitMod(dose, resp, model = "emax", addArgs = addArgs, bnds = bnds$emax);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
exponential = {fit <- DoseFinding::fitMod(dose, resp, model = "exponential", addArgs = addArgs, bnds = bnds$exponential);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
linlog = {fit <- DoseFinding::fitMod(dose, resp, model = "linlog", addArgs = addArgs, bnds = bnds$linlog);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
sigEmax = {fit <- DoseFinding::fitMod(dose, resp, model = "sigEmax", addArgs = addArgs, bnds = bnds$sigEmax);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
quadratic = {fit <- DoseFinding::fitMod(dose, resp, model = "quadratic", addArgs = addArgs, bnds = bnds$quadratic);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
betaMod = {fit <- DoseFinding::fitMod(dose, resp, model = "betaMod", addArgs = addArgs, bnds = bnds$betaMod);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
logistic = {fit <- DoseFinding::fitMod(dose, resp, model = "logistic", addArgs = addArgs, bnds = bnds$logistic);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
stop("The vector 'models' contains an incorrect entry")
)
};
for(j in 1:length(models)){result[j,1]<-round(length(resp)*log(minT0/RSS[j]),4)};
result[length(models)+1,1]<-max(result[1:length(models),1]);
return(result)
}
# SIMULATING the p-value
pLRcontrast<-function(dose, resp, models, off = 0.01 * max(dose), scal = 1.2 * max(dose), nsim = 1000, info = TRUE)
{
env1$data<-resp[order(dose)];
if(length(dose)!=length(resp)){stop("The length of the vectors 'dose' and 'resp' needs to be equal")};
env1$dose<-sort(unique(dose));
env1$weight<-c(); for(i in 1:length(env1$dose)){env1$weight[i]<-sum(env1$dose[i]==dose)/length(dose)};
env1$info=info;
if(min(env1$dose) < 0){stop("Only dose-levels >= 0 allowed")};
env1$off<-off;
if(env1$off <= 0){stop("'off' parameter needs to be positive")};
env1$scal<-scal;
if(env1$scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
if(nsim<=0 | nsim%%1!=0){stop("'nsim' needs to be a positive integer")};
env1$statistic<-sLRcontrast(dose=dose, resp=resp, models=models, off=env1$off, scal=env1$scal);
env1$Z<-matrix(data=rnorm(length(env1$dose)*nsim),nrow=length(env1$dose));
env1$Y<-matrix(data=rnorm(length(env1$dose)*nsim),nrow=length(env1$dose));
env1$L<-matrix(nrow=length(models),ncol=nsim);
env1$L_max<-c();
maxFkt<-list()
result<-matrix(nrow=length(models), ncol= 2); rownames(result)<-models; colnames(result)<-c("unadj-p", "adj-p");
for(j in 1:length(models))
{
switch(models[j],
linear={maxFkt[[j]]<-MAX_linear},
emax={maxFkt[[j]]<-MAX_emax},
exponential={maxFkt[[j]]<-MAX_exponential},
linlog={maxFkt[[j]]<-MAX_linlog},
sigEmax={maxFkt[[j]]<-MAX_sigEmax},
quadratic={maxFkt[[j]]<-MAX_quadratic},
betaMod={maxFkt[[j]]<-MAX_betaMod},
logistic={maxFkt[[j]]<-MAX_logistic},
stop("The vector 'models' contains an incorrect entry")
)
};
start<-Sys.time();
if(info==TRUE){cat(paste("% done \t", "nsim done\t", "secs passed\t", "secs remaining ~", collapse=""), sep="", fill=TRUE)};
for(j in 1:nsim){env1$i<-j; M1_max<-M1max(1);
for(k in 1:length(models)){env1$L[k,j]<-max(maxFkt[[k]](1)-M1_max,0)};
env1$L_max[j]<-max(env1$L[,j]);
if(info==TRUE){
passed <- round(difftime(Sys.time(),start,units="secs")[[1]],0);
percent<- round(j/nsim*100,2);
if(j%%10==0){env1$remaining<- round(passed/j*(nsim-j),0)};
cat(paste(rep.int("\b", 60),collapse=""),sep="");
cat(paste(percent, "\t\t", j, "\t\t", passed, "\t\t", env1$remaining, collapse=""), sep="");
flush.console()}
};
for(k in 1:length(models)){result[k,1]<-round(sum(env1$L[k,] >= env1$statistic[k])/nsim,4)};
for(k in 1:length(models)){result[k,2]<-round(sum(env1$L_max >= env1$statistic[k])/nsim,4)};
cat("\n");
return(result)}
Try the LRcontrast package in your browser
Any scripts or data that you put into this service are public.
LRcontrast documentation built on May 2, 2019, 5:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#SIMULATING QUANTILES:
#define variables:
#not changeable for the user:
env1<-new.env() #environment in which scal, off, i is saved
SCALE1<-c(1,1,1,1) #scaling of the parameters in fct optim (parscale)
LOW1<-c(0,0,0,0) #lower bounds for (phi1, phi2, theta1, theta2)
UP1<-c(pi,pi/2,6,6) #upper bounds for (phi1, phi2, theta) (LOW and UP needed in nlminb)
INTERVAL<-c(0, pi) #interval for phi (needed in optimize)
#starting values for nlminb in qLRcontrast and pLRcontrast:
#for the models linear, linlog, quadratic, emax, exponential:
initVal1<-matrix(nrow=3,ncol=3) #columns of the matrix are different starting values for (phi1, phi2, theta) in nlminb
initVal1[3,]<-rep(2,3)
initVal1[2,]<-rep(1,3)
initVal1[1,]<-c(0,pi/2,pi)
#for the models logistic and betaMod
initVal2<-matrix(nrow=4,ncol=8)
initVal2[4,]<-rep(2,8)
initVal2[3,]<-rep(2,8)
initVal2[2,]<-rep(1,8)
initVal2[1,]<-c(0,0.5,1,1.5,2,2.5,3,pi)
#for the model sigEmax
initVal3<-matrix(nrow=4,ncol=4)
initVal3[4,]<-rep(2,4)
initVal3[3,]<-rep(2,4)
initVal3[2,]<-rep(1,4)
initVal3[1,]<-c(0,1,2,pi)
#define the stochastic process:
#factor (-1), because nlminb minimizes functions instead of maximizing ( max{f(x)} = -min{-f(x)} )):
#parametrization with polar coordinates:
#here: beta0 = cos(phi), sigma = sin(phi)
WS1<-function(phi){return(sum(sqrt(env1$weight)*(cos(phi)*env1$Z[,env1$i]+sqrt(2)*sin(phi)*env1$Y[,env1$i]))/sqrt(sum(env1$weight*(cos(phi)^2+2*sin(phi)^2))))}
#here: beta0 = cos(phi1)*cos(phi2), beta1 = sin(phi2), sigma = sin(phi1)*cos(phi2)
WS2_linear<-function(x){return(.C("WS2linear", phi1=as.double(x[1]), phi2=as.double(x[2]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_emax<-function(x){return(.C("WS2emax", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_exponential<-function(x){return(.C("WS2exponential", phi1=as.double(x[1]), phi2=as.double(x[2]), delta=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_linlog<-function(x){return(.C("WS2linlog", phi1=as.double(x[1]), phi2=as.double(x[2]), off=as.double(env1$off), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_sigEmax<-function(x){return(.C("WS2sigEmax", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), h=as.double(x[4]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_quadratic<-function(x){return(.C("WS2quadratic", phi1=as.double(x[1]), phi2=as.double(x[2]), b=as.double(x[3]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_betaMod<-function(x){return(.C("WS2betaMod", phi1=as.double(x[1]), phi2=as.double(x[2]), delta1=as.double(x[3]), delta2=as.double(x[4]), scal=as.double(env1$scal), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
WS2_logistic<-function(x){return(.C("WS2logistic", phi1=as.double(x[1]), phi2=as.double(x[2]), ed50=as.double(x[3]), delta=as.double(x[4]), weight=as.double(env1$weight), dose=as.double(env1$dose), Z=as.double(env1$Z[,env1$i]), Y=as.double(env1$Y[,env1$i]), doselength = as.integer(length(env1$dose)))$phi1)}
#maximizing the functions:
#maximize WS1 over INTERVAL=[0,2*pi] (M1) (equal for all models)
M1max<-function(x){return(max(optimize(WS1, interval=INTERVAL, maximum=TRUE)$objective, 0)^2)}
#maximize WS2 for the different models over M2. Maximize with nlminb with different starting values and take the maximum.
MAX_linear<-function(x){max(max(0,(-1)*nlminb(initVal1[1:2,1], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,2], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,3], WS2_linear, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_emax<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_emax, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_exponential<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_exponential, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_linlog<-function(x){max(max(0,(-1)*nlminb(initVal1[1:2,1], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,2], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[1:2,3], WS2_linlog, lower=LOW1[1:2], upper=UP1[1:2], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_quadratic<-function(x){max(max(0,(-1)*nlminb(initVal1[,1], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,2], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal1[,3], WS2_quadratic, lower=LOW1[1:3], upper=UP1[1:3], control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_betaMod<-function(x){max(max(0,(-1)*nlminb(initVal2[,1], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,2], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,3], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,4], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,5], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,6], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,7], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,8], WS2_betaMod, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_logistic<-function(x){max(max(0,(-1)*nlminb(initVal2[,1], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,2], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,3], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,4], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,5], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,6], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,7], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal2[,8], WS2_logistic, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
MAX_sigEmax<-function(x){max(max(0,(-1)*nlminb(initVal3[,1], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,2], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,3], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2,
max(0,(-1)*nlminb(initVal3[,4], WS2_sigEmax, lower=LOW1, upper=UP1, control=list(iter.max=5000,eval.max=5000))$objective)^2)}
#with that function the user can simulate the quantiles:
qLRcontrast<-function(dose, probs, models, weight = rep(1/length(dose), length(dose)), off = 0.01 * max(dose), scal = 1.2 * max(dose), nsim = 10000, info = TRUE)
{
env1$weight<-weight;
if(sum(env1$weight)!=1){stop("The sum of 'weight' needs to be 1")};
if(min(env1$weight)<=0){stop("Only positive values in 'weight' allowed")};
env1$dose<-dose;
if(min(env1$dose) < 0){stop("Only dose-levels >= 0 allowed")};
if(length(env1$dose)!=length(env1$weight)){stop("Lengths of the vectors 'dose' and 'weight' need to be equal")};
env1$probs<-sort(probs);
if(env1$probs[1]<0 | 1<max(env1$probs)){stop("'probs' needs to be a numeric vector of probabilities with values in [0,1]")};
env1$off<-off;
if(env1$off <= 0){stop("'off' parameter needs to be positive")};
env1$scal<-scal;
if(env1$scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
if(nsim<=0 | nsim%%1!=0){stop("'nsim' needs to be a positive integer")};
env1$Z<-matrix(data=rnorm(length(dose)*nsim),nrow=length(dose));
env1$Y<-matrix(data=rnorm(length(dose)*nsim),nrow=length(dose));
env1$L<-matrix(nrow=length(models),ncol=nsim);
L_max<-c();
QUANTILE<-matrix(nrow=length(models)+1, ncol=length(env1$probs)); colnames(QUANTILE)<-env1$probs; rownames(QUANTILE)<-c(models, "max");
maxFkt<-list()
for(j in 1:length(models))
{
switch(models[j],
linear={maxFkt[[j]]<-MAX_linear},
emax={maxFkt[[j]]<-MAX_emax},
exponential={maxFkt[[j]]<-MAX_exponential},
linlog={maxFkt[[j]]<-MAX_linlog},
sigEmax={maxFkt[[j]]<-MAX_sigEmax},
quadratic={maxFkt[[j]]<-MAX_quadratic},
betaMod={maxFkt[[j]]<-MAX_betaMod},
logistic={maxFkt[[j]]<-MAX_logistic},
stop("The vector 'models' contains an incorrect entry")
)
};
start<-Sys.time();
if(info==TRUE){cat(paste("% done \t", "nsim done\t", "secs passed\t", "secs remaining ~", collapse=""), sep="", fill=TRUE)};
for(j in 1:nsim){env1$i<-j; M1_max<-M1max(1);
for(k in 1:length(models)){env1$L[k,j]<-maxFkt[[k]](1)-M1_max};
L_max[j]<-max(env1$L[,j])
if(info==TRUE){
passed <- round(difftime(Sys.time(),start,units="secs")[[1]],0);
percent<- round(j/nsim*100,2);
if(j%%10==0){env1$remaining<- round(passed/j*(nsim-j),0)};
cat(paste(rep.int("\b", 60),collapse=""),sep="");
cat(paste(percent, "\t\t", j, "\t\t", passed, "\t\t", env1$remaining, collapse=""), sep="");
flush.console()}
};
for(j in 1:length(models)){QUANTILE[j,]<-quantile(env1$L[j,],env1$probs)};
QUANTILE[length(models)+1,]<-quantile(L_max,env1$probs);
cat("\n");
return(QUANTILE)}
#test statistic (2.9):
sLRcontrast<-function(dose, resp, models, off = 0.01 * max(dose), scal= 1.2 * max(dose))
{
if(length(dose)!=length(resp)){stop("The length of the vectors 'dose' and 'resp' needs to be equal")};
if(min(dose) < 0){stop("Only dose-levels >= 0 allowed")};
if(off <= 0){stop("'off' parameter needs to be positive")};
if(scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
result <- matrix(nrow=length(models)+1, ncol= 1); rownames(result)<-c(models, "max"); colnames(result)<-c("statistic");
RSS <- c();
addArgs <- list(); addArgs$off <- off; addArgs$scal <- scal;
bnds <- DoseFinding::defBnds(mD = max(dose));
minT0<-sum((resp-1/length(resp)*sum(resp))^2);
for(j in 1:length(models))
{
switch(models[j],
linear = {fit <- DoseFinding::fitMod(dose, resp, model = "linear", addArgs = addArgs, bnds = bnds$linear);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
emax = {fit <- DoseFinding::fitMod(dose, resp, model = "emax", addArgs = addArgs, bnds = bnds$emax);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
exponential = {fit <- DoseFinding::fitMod(dose, resp, model = "exponential", addArgs = addArgs, bnds = bnds$exponential);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
linlog = {fit <- DoseFinding::fitMod(dose, resp, model = "linlog", addArgs = addArgs, bnds = bnds$linlog);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
sigEmax = {fit <- DoseFinding::fitMod(dose, resp, model = "sigEmax", addArgs = addArgs, bnds = bnds$sigEmax);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
quadratic = {fit <- DoseFinding::fitMod(dose, resp, model = "quadratic", addArgs = addArgs, bnds = bnds$quadratic);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
betaMod = {fit <- DoseFinding::fitMod(dose, resp, model = "betaMod", addArgs = addArgs, bnds = bnds$betaMod);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
logistic = {fit <- DoseFinding::fitMod(dose, resp, model = "logistic", addArgs = addArgs, bnds = bnds$logistic);
if(fit$coefs[[2]] > 0){RSS[j] <- fit$RSS} else{RSS[j] <- minT0}},
stop("The vector 'models' contains an incorrect entry")
)
};
for(j in 1:length(models)){result[j,1]<-round(length(resp)*log(minT0/RSS[j]),4)};
result[length(models)+1,1]<-max(result[1:length(models),1]);
return(result)
}
# SIMULATING the p-value
pLRcontrast<-function(dose, resp, models, off = 0.01 * max(dose), scal = 1.2 * max(dose), nsim = 1000, info = TRUE)
{
env1$data<-resp[order(dose)];
if(length(dose)!=length(resp)){stop("The length of the vectors 'dose' and 'resp' needs to be equal")};
env1$dose<-sort(unique(dose));
env1$weight<-c(); for(i in 1:length(env1$dose)){env1$weight[i]<-sum(env1$dose[i]==dose)/length(dose)};
env1$info=info;
if(min(env1$dose) < 0){stop("Only dose-levels >= 0 allowed")};
env1$off<-off;
if(env1$off <= 0){stop("'off' parameter needs to be positive")};
env1$scal<-scal;
if(env1$scal < max(dose)){stop("'scal' parameter needs to be >= max(dose)")};
if(length(models)!=length(unique(models))){stop("Only one list entry allowed for each model")};
if(nsim<=0 | nsim%%1!=0){stop("'nsim' needs to be a positive integer")};
env1$statistic<-sLRcontrast(dose=dose, resp=resp, models=models, off=env1$off, scal=env1$scal);
env1$Z<-matrix(data=rnorm(length(env1$dose)*nsim),nrow=length(env1$dose));
env1$Y<-matrix(data=rnorm(length(env1$dose)*nsim),nrow=length(env1$dose));
env1$L<-matrix(nrow=length(models),ncol=nsim);
env1$L_max<-c();
maxFkt<-list()
result<-matrix(nrow=length(models), ncol= 2); rownames(result)<-models; colnames(result)<-c("unadj-p", "adj-p");
for(j in 1:length(models))
{
switch(models[j],
linear={maxFkt[[j]]<-MAX_linear},
emax={maxFkt[[j]]<-MAX_emax},
exponential={maxFkt[[j]]<-MAX_exponential},
linlog={maxFkt[[j]]<-MAX_linlog},
sigEmax={maxFkt[[j]]<-MAX_sigEmax},
quadratic={maxFkt[[j]]<-MAX_quadratic},
betaMod={maxFkt[[j]]<-MAX_betaMod},
logistic={maxFkt[[j]]<-MAX_logistic},
stop("The vector 'models' contains an incorrect entry")
)
};
start<-Sys.time();
if(info==TRUE){cat(paste("% done \t", "nsim done\t", "secs passed\t", "secs remaining ~", collapse=""), sep="", fill=TRUE)};
for(j in 1:nsim){env1$i<-j; M1_max<-M1max(1);
for(k in 1:length(models)){env1$L[k,j]<-max(maxFkt[[k]](1)-M1_max,0)};
env1$L_max[j]<-max(env1$L[,j]);
if(info==TRUE){
passed <- round(difftime(Sys.time(),start,units="secs")[[1]],0);
percent<- round(j/nsim*100,2);
if(j%%10==0){env1$remaining<- round(passed/j*(nsim-j),0)};
cat(paste(rep.int("\b", 60),collapse=""),sep="");
cat(paste(percent, "\t\t", j, "\t\t", passed, "\t\t", env1$remaining, collapse=""), sep="");
flush.console()}
};
for(k in 1:length(models)){result[k,1]<-round(sum(env1$L[k,] >= env1$statistic[k])/nsim,4)};
for(k in 1:length(models)){result[k,2]<-round(sum(env1$L_max >= env1$statistic[k])/nsim,4)};
cat("\n");
return(result)}
Try the LRcontrast package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.