diff --git a/DESCRIPTION b/DESCRIPTION index d12ec83..bcf4ebd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: ToxicR Type: Package Title: Analyzing Toxicology Dose-Response Data -Version: 22.8.1.0.4 -Date: 2022-10-25 +Version: 22.8.1.0.5 +Date: 2022-11-10 Authors@R: c( person(given = "Matt", diff --git a/R/MAdensity_plot.R b/R/MAdensity_plot.R index 8208f39..5db2f9a 100644 --- a/R/MAdensity_plot.R +++ b/R/MAdensity_plot.R @@ -82,11 +82,7 @@ MAdensity_plot <- function (A){ - # Q <- t(Q) - # - # me <- colMeans(Q) - # lq <- apply(Q,2,quantile, probs = qprob) - # uq <- apply(Q,2,quantile, probs = 1-qprob) + Dens = density(temp,cut=c(max(doses))) # what is this 0.4 means? Scale? @@ -128,37 +124,7 @@ MAdensity_plot <- function (A){ idx <- sample(1:length(fit_idx), length(A$Individual_Model_1$mcmc_result$BMD_samples),replace=TRUE,prob=A$posterior_probs) df<-NA - # - # m1<-A$Individual_Model_1 - # c1<-m1$mcmc_result$BMD_samples - # - # - # m2<-A$Individual_Model_2 - # c2<-m2$mcmc_result$BMD_samples - # - # m3<-A$Individual_Model_3 - # c3<-m3$mcmc_result$BMD_samples - # - # m4<-A$Individual_Model_4 - # c4<-m4$mcmc_result$BMD_samples - # - # m5<-A$Individual_Model_5 - # c5<-m5$mcmc_result$BMD_samples - # - # m6<-A$Individual_Model_6 - # c6<-m6$mcmc_result$BMD_samples - # - # m7<-A$Individual_Model_7 - # c7<-m7$mcmc_result$BMD_samples - # - # m8<-A$Individual_Model_8 - # c8<-m5$mcmc_result$BMD_samples - # - # m9<-A$Individual_Model_9 - # c9<-m9$mcmc_result$BMD_samples - # - # combine_samples<-data.frame(cbind(c1,c2,c3,c4,c5,c6,c7,c8,c9)) - + # How should I initialize this? for (i in 1:length(fit_idx)){ @@ -196,37 +162,7 @@ MAdensity_plot <- function (A){ t_combine3<-t_combine2 %>% filter(as.numeric(X3)>0.05 , as.numeric(X1)!=Inf) - - - - # # Update to change color - # - # step1<- p<-ggplot()+ - # stat_density_ridges(data=t_combine3, aes(x=X1, y=fct_reorder(X2,X3,.desc=T),group=X2, - # fill=cut(X3,c(0,0.999,1))), - # calc_ecdf=TRUE,quantiles=c(0.025,0.975),na.rm=T,quantile_lines=T,scale=0.9)+ - # scale_fill_manual(name="X3", values=c("(0,0.999]"="grey", "(0.999,1]"="blue")) - # - # - # step2<-step1+ - # xlim(c(0,quantile(t_combine$X1,0.999)))+ - # geom_vline(xintercept = A$bmd[1],linetype="longdash")+ - # labs(x="Dose Level (Dotted Line : MA BMD)",y="", title="Density plots for each fitted model (Fit type: MCMC)")+theme_classic() - # - # step2 - # - # step3<-step2+ - # stat_density_ridges(data=t_combine3, aes(x=X1, y=fct_reorder(X2,X3,.desc=T), group=X2, fill = factor(stat(quantile))), - # geom = "density_ridges_gradient", - # calc_ecdf = TRUE, quantiles = c(0.025, 0.975), scale=0.9)+ - # scale_fill_manual( - # name = "Probability", - # values = c("NA","NA","NA", "NA", "NA") - # ) - # - # - # step3 - # + p<-ggplot()+ @@ -241,15 +177,6 @@ MAdensity_plot <- function (A){ p2<-p+theme(legend.position="none", axis.text.y=element_text(size=12)) - # stat_density_ridges(data=t_combine3, aes(x=X1, y=fct_reorder(X2,X3,.desc=T), group=X2,fill = factor(stat(quantile))), - # geom = "density_ridges_gradient", - # calc_ecdf = TRUE, quantiles = c(0.025, 0.975),scale=0.9 - # )+ scale_fill_manual( - # name = "factor(stat(quantile))", - # values = c("#FF0000A0", "NA", "#FF0000A0"), - # labels = c("(0, 0.025]", "(0.025, 0.975]", "(0.975, 1]") - # )+ - p2 return(p2) # Return output @@ -258,8 +185,6 @@ MAdensity_plot <- function (A){ } -# Laplace case- Is it even possible to do this? -# No we don't need this part .plot.density.BMDdichotomous_MA_maximized<-function(A){ t_1 <- t_2 <- t_3 <- t_4 <- t_5 <- t_6 <- t_7 <- t_8 <- t_9 <- c3 <- X1 <- X2 <- X3 <- NULL @@ -406,152 +331,56 @@ MAdensity_plot <- function (A){ for (i in fit_idx){ # Loop for the model fit<-A[[i]] - test_doses <- seq(min(doses),max(doses)*1.03,(max(doses)*1.03-min(doses))/100) - #probs <- (0.5+fit$data[,2,drop=T])/(1.0 + fit$data[,3,drop=T]) - + + temp <- fit$mcmc_result$BMD_samples[!is.nan(fit$mcmc_result$BMD_samples)] + temp <- temp[!is.infinite(temp)] - if (fit$model=="hill"){ - Q <- apply(fit$mcmc_result$PARM_samples,1,.cont_hill_f, d=test_doses) - - } - if (fit$model=="exp-3"){ - Q <- apply(fit$mcmc_result$PARM_samples,1,.cont_exp_3_f, d=test_doses) - - } - if (fit$model=="exp-5"){ - Q <- apply(fit$mcmc_result$PARM_samples,1,.cont_exp_5_f, d=test_doses) + # Try to run the density function. + # If there is a problem, return null + + temp_density<-data.frame(matrix(0,length(temp),3)) + temp_density[,2]=substr(fit$full_model,8,999) + temp_density[,1]=temp + temp_density[,3]=A$posterior_probs[i] - } - if (fit$model=="power"){ - Q <- apply(fit$mcmc_result$PARM_samples,1,.cont_power_f, d=test_doses) + # assign(paste("t",i,sep="_"),temp_density) - } - if (fit$model=="FUNL"){ - Q <- apply(fit$mcmc_result$PARM_samples,1,.cont_FUNL_f, d=test_doses) + t<-temp_density - } - - - temp <- fit$mcmc_result$BMD_samples[!is.nan(fit$mcmc_result$BMD_samples)] - temp <- temp[!is.infinite(temp)] - - - - # Q <- t(Q) - # - # me <- colMeans(Q) - # lq <- apply(Q,2,quantile, probs = qprob) - # uq <- apply(Q,2,quantile, probs = 1-qprob) - - Dens = density(temp,cut=c(max(doses))) - # what is this 0.4 means? Scale? - - # normalize it?-- We don't need it actually here - # Dens$y = Dens$y/max(Dens$y) * max(probs) - # temp = which(Dens$x < max(doses)) - # D1_y = Dens$y[temp] - # D1_x = Dens$x[temp] - - - - # Do I need to stack up the dataset? - - - temp_density<-data.frame(matrix(0,length(temp),3)) - temp_density[,2]=substr(fit$full_model,8,999) - temp_density[,1]=temp - temp_density[,3]=A$posterior_probs[i] - - # assign(paste("t",i,sep="_"),temp_density) - - t<-temp_density - - t_combine<-rbind(t_combine,t) - + t_combine<-rbind(t_combine,t) + } t_combine<-t_combine[-1,] + idx <- sample(1:length(fit_idx), length(A$Individual_Model_1$mcmc_result$BMD_samples), + replace=TRUE,prob=A$posterior_probs) - - # This part is needed to get MA density plots - # Model ~xxx.. - # Number of model should be fixed - - # Grep function -- index - fit_idx - - - - idx <- sample(1:length(fit_idx), length(A$Individual_Model_1$mcmc_result$BMD_samples),replace=TRUE,prob=A$posterior_probs) - - - - - # rbind(A[[1]]$mcmc_result$BMD_samples) - df<-NA - - - # How should I initialize this? + ## for (i in 1:length(fit_idx)){ m<-A[[i]] c<-m$mcmc_result$BMD_samples - df<-data.frame(cbind(df,c)) + df<-cbind(df,c) } - # - # m1<-A$Individual_Model_1 - # c1<-m1$mcmc_result$BMD_samples - # - # - # m2<-A$Individual_Model_2 - # c2<-m2$mcmc_result$BMD_samples - # - # m3<-A$Individual_Model_3 - # c3<-m3$mcmc_result$BMD_samples - # - # m4<-A$Individual_Model_4 - # c4<-m4$mcmc_result$BMD_samples - # - # m5<-A$Individual_Model_5 - # c5<-m5$mcmc_result$BMD_samples - - # This part should be dynamically assigned... - # combine_samples2<-data.frame(cbind(c1,c2,c3,c4,c5)) - - df_samples<-data.frame(df[,-1]) - - - - - - # Select MA values - - BMD_MA<-matrix(NA,length(A$Individual_Model_1$mcmc_result$BMD_samples),1) - - for (i in 1:length(A$Individual_Model_1$mcmc_result$BMD_samples)){ - j<-sample(nrow(df_samples), size=1, replace=TRUE) - BMD_MA[i,1]<-df_samples[j,idx[i]] + # Compute the model average density + df <- as.matrix(df[,-1]) + result_idx = sample(1:ncol(df),dim(df)[1],replace=T,prob= A$posterior_probs ) + BMD_MA = matrix(NA,dim(df)[1],3) + for (ii in 1:(dim(BMD_MA)[1])){ + BMD_MA[ii,1] = df[ii,result_idx[ii]] } - - - BMD_MA<-data.frame(BMD_MA) - - t_ma<-BMD_MA %>% mutate(X1=BMD_MA,X2="Model Average",X3=1) - #BMD_CDF should be shown here - it - t_ma<-t_ma %>% select(X1,X2,X3) - - - - t_combine2<-rbind(t_combine,t_ma) - - - t_combine3<-t_combine2 %>% filter(as.numeric(X3)>0.05 , as.numeric(X1)!=Inf) - - - # John's comment- I want to fill the color as + BMD_MA[,2] = "Model Average" + BMD_MA[,3] = 1 + #clean up t_combine for the plot + t_combine <- t_combine %>% filter(as.numeric(X3)>0.05 , as.numeric(X1)!=Inf) + t_combine <-rbind(t_combine,BMD_MA) + t_combine3 <- t_combine %>% filter(!is.infinite(as.numeric(X3)), + !is.na(as.numeric(X1))) + + #make plot p<-ggplot()+ stat_density_ridges(data=t_combine3, aes(x=X1, y=fct_reorder(X2,X3,.desc=T),group=X2 ,alpha=sqrt(X3), fill=cut(X3,c(0,0.99,1))), calc_ecdf=TRUE,quantiles=c(0.025,0.975),na.rm=T,quantile_lines=T,scale=0.9)+ @@ -568,7 +397,6 @@ MAdensity_plot <- function (A){ return(p2) # Return output - } diff --git a/R/cleveland_plot.R b/R/cleveland_plot.R index 0c6b53c..4e5681c 100644 --- a/R/cleveland_plot.R +++ b/R/cleveland_plot.R @@ -108,7 +108,7 @@ cleveland_plot <- function (A){ bmd_ind_df2<-data.frame(bmd_ind_df[which(!is.na(bmd_ind_df[,1])),]) out<-ggplot()+ - geom_point(data=bmd_ind_df2, aes(x=as.numeric(X1), y=fct_reorder(X4,as.numeric(X5),.desc=T),size=(sqrt(as.numeric(X5)))), color="red")+ + geom_point(data=bmd_ind_df2, aes(x=as.numeric(X1), y=fct_reorder(X4,as.numeric(X5),.desc=T),size=(sqrt(as.numeric(X5)+0.01))), color="red")+ #scale_colour_gradient(low = "gray", high = "black")+ geom_vline(xintercept=as.numeric(bmd_ind_df2$X1[length(fit_idx)+1]), linetype="dotted")+ theme_minimal()+ @@ -165,9 +165,9 @@ cleveland_plot <- function (A){ bmd_ind_df2<-data.frame(bmd_ind_df[which(!is.na(bmd_ind_df[,1])),]) out<-ggplot()+ - geom_point(data=bmd_ind_df2, aes(x=as.numeric(X1), y=fct_reorder(X4,as.numeric(X5),.desc=T),size=(sqrt(as.numeric(X5)))), color="red")+ + geom_point(data=bmd_ind_df2, aes(x=as.numeric(X1), y=fct_reorder(X4,as.numeric(X5),.desc=T),size=(sqrt(0.01+as.numeric(X5)))), color="red")+ #scale_colour_gradient(low = "gray", high = "black")+ - geom_vline(xintercept=as.numeric(bmd_ind_df2$X1[length(fit_idx)+1]), linetype="dotted")+ + #geom_vline(xintercept=as.numeric(bmd_ind_df2$X1[length(fit_idx)+1]), linetype="dotted")+ theme_minimal()+ labs(x="Dose Level",y="", title="BMD Estimates by Each Model (Sorted by Posterior Probability)",size="Posterior Probability")+ theme(legend.position="none")+ diff --git a/vignettes/Continuous.Rmd b/vignettes/Continuous.Rmd index e699ab7..c537e44 100755 --- a/vignettes/Continuous.Rmd +++ b/vignettes/Continuous.Rmd @@ -26,7 +26,7 @@ This file shows ToxicR for dichotomous benchmark dose analyses using both single cont_data <- matrix(0,nrow=5,ncol=4) colnames(cont_data) <- c("Dose","Mean","N","SD") cont_data[,1] <- c(0,50,100,200,400) -cont_data[,2] <- c(5.26,5.76,6.13,8.24,9.23) +cont_data[,2] <- c(5.26,5.76,7.13,9.24,9.23) cont_data[,3] <- c(20,20,20,20,20) cont_data[,4]<- c(2.23,1.47,2.47,2.24,1.56) Y <- cont_data[,2:4] @@ -243,28 +243,6 @@ plot(poly_sd) ## Model Averaging -Beyond single model fits, one can also use model averging for continuous models. -As a default 10 models are fit. Here, the 'hill' and 'power' options are fit under -the 'normal' and 'normal-ncv' distribution assumptions. Additionally to the -'normal' and 'normal-ncv' distributions the 'exp-3' and 'exp-5' also use the -'lognormal' distribution. Note that the polynomial model is not used because -of constraint issues. - -```{r run_laplace_MA, fig.height = 5, fig.width = 6} -ma_sd_mcmc <- ma_continuous_fit(cont_data[,"Dose"],Y, fit_type="mcmc", - BMD_TYPE="sd",BMR = 0.5,samples = 50000) - -ma_sd_laplace <- ma_continuous_fit(cont_data[,"Dose"],Y, fit_type="laplace", - BMD_TYPE="sd",BMR = 0.5) - -plot(ma_sd_mcmc) -cleveland_plot(ma_sd_laplace) -MAdensity_plot(ma_sd_mcmc) -``` - -Here, 'cleveland_plot' and 'MAdensity_plot' can be used in the same way as in -dichotomous plot. Like the dichotomous case, we can create a prior list for our model -average. ```{r run_laplace_MA_2} prior <- create_prior_list(normprior(0,1,-100,100),