In cdv04/ACTR:

---

knitr::opts_chunk$set(echo = TRUE, message=FALSE)
#set.seed(1234) # in order to obtain same results
# load package

library(dplyr)
library(ACTpckg)
sessionInfo()

randNum <- 1234

B <- 1000 # number of bootstrap samples

coef <- 0 # coef for detecting outliers

nbED <- 6

ptm <- proc.time()


# load data
data(cipr)

# check step
check_all(cipr) # print levels of different variable to check orthography into a txt file
nb_ED(cipr, Class, nbED) # print levels of Class having at least 6 EDR10 and 6ED50


# keep data of Classes having at least 6 EDR10 and 6 ED50
ciprKP0 <- subset_data(cipr, Class,nbED)

# write cipr_kp data into xlsx file for SSWD check
# library(xlsx)
# write.xlsx(cipr_kp,"D:\\ProjetsIRSN\\ACT\\ACT_v2\\VerifParaSSWD\\cipr_kp.xlsx" )

ciprKP0.Ac <- filter(ciprKP0, DoseType=="Acute")
nb_outliers(ciprKP0.Ac,Class,coef) 
show_outliers(ciprKP0.Ac,Class,coef) 


ciprKP0.Ch <- filter(ciprKP0, DoseType=="Chronic")
nb_outliers(ciprKP0.Ch,Class,coef) 
show_outliers(ciprKP0.Ch,Class,coef) 

ciprKP1.Ac<-rm_outliers(ciprKP0.Ac,coef)
ciprKP1.Ch<-rm_outliers(ciprKP0.Ch,coef)

ciprKP1 <- rbind(ciprKP1.Ac,ciprKP1.Ch)

nb_ED(ciprKP1, Class, nbED) # print levels of Class having at least 6 EDR10 and 6ED50 in a txt file


# keep data of Classes having at least 6 EDR10 and 6 ED50

ciprKP <- subset_data(ciprKP1, Class,nbED)

set.seed(randNum ) # in order to obtain same results


###############################################################################
##  Selecting best model for prediction of mu_C and Sigma_C, by bootstrap and
##  cross validation
##


# bootstrap sampling of the cipr_kp data
ciprKP.bt <- boot_data(ciprKP,Class,B)


# estimation of weighted (on species) log10 Mu_A, Mu_C, Sigma_A and Sigma_C of
#bootstrap sample
PAC_bt <- lapply(ciprKP.bt,est_PAC,Class)



# estimation of prediction error of several models for each of the bootstrap
#sample
errMod_mu <- lapply(PAC_bt,msep_mu)
errMod_sigma <- lapply(PAC_bt,msep_sigma)


# choose the best model (most little Mean Square Prediction Error) by Cross
#Validation
bm_mu <- bestModel_mu(errMod_mu)
bm_mu
bm_sigma <- bestModel_sigma(errMod_sigma)
bm_sigma

################################################################################
##  Generating EDR10_ACT
##


# estimation of observed weighted parameters
PAC_obs <- est_PAC(ciprKP,Class)
PAC_obs
# predicted Chronic parameters (using best models and observed acute parameters)
PC_pred <- predict_PC(PAC_obs,bm_mu, bm_sigma)
PC_pred
coef_Mod <- coef_reg (PAC_obs,bm_mu, bm_sigma)
coef_Mod
sd_res <- sd_res_reg (PAC_obs,bm_mu, bm_sigma)
sd_res
# generating edr10_act
ciprKP.ac <- filter(ciprKP, DoseType=="Acute")
PA_obs<-arrange(PAC_obs[,1:3],PAC_obs[,1])
PC_pred <- arrange(PC_pred,PC_pred[,1])
PAC_gen <- cbind(PA_obs,PC_pred[,-1] )

edr10_act <- gen_edr10(PAC_gen,ciprKP.ac)

###############################################################################
##  plotting Empirical Cumulative Distribution of edr10 observed
##  and edr10 generated by act
##

# spliting edr10_obs from origin data
edr10_obs <- split_edr10_obs(ciprKP, Class)

# calculating ploting position
edr10_obs <- plotting_position (edr10_obs)
edr10_act <- plotting_position (edr10_act)


# plot
edr10_all <- rbind(edr10_obs, edr10_act)

plot_ecdf <- plotECDF(edr10_all)
plot_ecdf

##############################################################################
##  plotting Log Normal  Cumulative Distribution of edr10 observed
##  and edr10 generated by act
##


plot_ln <- plotLN(edr10_all)
plot_ln

# in order to obtain same results

###############################################################################
##  Estimating HDR5 obs and IC from empirical distribution
## NB : EO stands for "Empirical &  Observed"
#hdr5_LNO <- Hdr_LN(edr10_obs,0.05)
set.seed(randNum )

hdr5_LNO <-IC_HdrLN(edr10_obs, 0.05,B)
hdr5_LNO

#hdr5_LNO2 <-IC_HdrLN2(edr10_obs, 0.05,1000) # test IC sans pondération Mu et Sigma comme ds SSWD

###############################################################################
##  Estimating HDR5 ACT (from EDR10_ACT) from Log Normal distribution
##
set.seed(randNum )

hdr5_LNA <- IC_HdrLN_act(ciprKP, Class, bm_mu, bm_sigma, 0.05,B)
hdr5_LNA

###############################################################################
##  Comparing hdr5_obs and hdr5_act
##

myhdr5<-resume_hdr(hdr5_LNO,hdr5_LNA)
myhdr5

plot_hdr5 <- view_hdr(hdr5_LNO, hdr5_LNA) # IC avec les moy et sd pondérés
plot_hdr5

 # in order to obtain same results

###############################################################################
##  Estimating HDR50 obs and IC from empirical distribution
## NB : EO stands for "Empirical &  Observed"
#hdr50_LNO <- Hdr_LN(edr10_obs,0.5)
set.seed(randNum )
hdr50_LNO <- IC_HdrLN(edr10_obs, 0.5,B)
hdr50_LNO
#hdr50_LNO2 <- IC_HdrLN2(edr10_obs, 0.5,1000) # test IC sans pondération Mu et Sigma comme ds SSWD

###############################################################################
##  Estimating HDR5 ACT (from EDR10_ACT) from log normal distribution
##
set.seed(randNum )
hdr50_LNA <- IC_HdrLN_act(ciprKP, Class, bm_mu, bm_sigma, 0.5,B)
hdr50_LNA 

###############################################################################
##  Comparing hdr5_obs and hdr5_act
##

myhdr50<-resume_hdr(hdr50_LNO,hdr50_LNA)
myhdr50

plot_hdr50 <-view_hdr(hdr50_LNO, hdr50_LNA) # IC avec les moy et sd pondérés
plot_hdr50

# plot_hdr50_bis <-view_hdr(hdr50_EO2, hdr50_EA) # IC sans les moy et sd pondérés
# plot_hdr50_bis

proc.time() - ptm

cdv04/ACTR documentation built on May 13, 2019, 2:42 p.m.