---
title: "Main figures for article on 20 years of the npde adventure, for AAPS J"
author: "Emmanuelle Comets"
date: "19/02/2021"
output:
  pdf_document: default
  html_notebook: default
  html_document:
    df_print: paged
---

## Introduction

This notebook is an R Markdown document containing R code for key figures in the paper **Developing tools to evaluate non-linear mixed effect models: 20 years on the npde adventure**, by Emmanuelle Comets and France Mentré, submitted to AAPS Journal for the special edition *Women in Science*. It is bundled in a Zenodo archive along with all the files necessary to reproduce the code and figures herein.

This document can be used to reproduce the main figures in the paper. It includes the citation data used in the article, and shows the data management applied to extract relevant keywords. A number of libraries are required for this code to run (see first code chunk).

At the end, we also show how to compute the ndpe for the viral load example used in the paper, producing the default graphs from the new npde 3.0 package. The npde package submitted to CRAN on the same day is included in the Zenodo archive (Linux version).

## Loading libraries

```{r setup}
knitr::opts_chunk$set(echo = TRUE)
# Libraries
library(gridExtra)
library(ggplot2)
library(grid)
library(RColorBrewer)
library(dplyr)
library(devtools)
library(mclust)

# Reading Excel files (?)
library("readxl")

# Word count libraries
library("tm")
library("SnowballC")
library("wordcloud")

```

## Most cited articles 

### Database of citations for 5 most cited articles

The citation data was downloaded from Web of Science on October 1st, 2020, for the 5 foundation articles mentioned in the submitted article.

```{r citationData}
cite.list<-dir(pattern="citations")
cite.data<-vector(length=length(cite.list),mode="list")
for(ipub in 1:length(cite.list)) {
  cite.data[[ipub]]<- read_excel(cite.list[ipub])
}
nampub<-unlist(strsplit(cite.list,"citations_"))
nampub<-unlist(strsplit(nampub,".xls"))
nampub<-nampub[nampub!=""]
```

Checking autocitations for each article

```{r autoCitation}
for(ipub in 1:length(cite.list)) {
  datcite<-cite.data[[ipub]]
  id1<-grep("Comets",datcite$Authors)
  id2<-grep("Mentr",datcite$Authors)
  id3<-grep("Brendel",datcite$Authors)
  id<-unique(c(id1,id2,id3))
  cat("Publication:",nampub[ipub],"\n")  
  cat("Number of citations:",dim(datcite)[1],"including",length(id),"auto-citations \n")
  cat("        number of articles:",dim(datcite[datcite$`Document Type`=="Article",])[1],"\n")
  print(table(datcite$`Publication Type`))
  print(table(datcite$`Document Type`))
}
```

### Data selection

- Tabulating data type and citations
- create dataset with all articles and letters

```{r numberCitations}
# Create dataset with all articles and letters, for all articles for the moment
fulldat<-books<-NULL
for(ipub in 1:length(cite.list)) {
  datcite<-cite.data[[ipub]]
  id1<-grep("Comets",datcite$Authors)
  id2<-grep("Mentr",datcite$Authors)
#  id3<-grep("Brendel",datcite$Authors)
  id<-unique(c(id1,id2))
  cat("Publication:",nampub[ipub],"\n")  
  cat("Number of citations:",dim(datcite)[1],"including",length(id),"auto-citations \n")
  cat("        number of articles:",dim(datcite[datcite$`Document Type`=="Article",])[1],"\n")
  print(table(datcite$`Publication Type`))
  print(table(datcite$`Document Type`))
  dat1<-datcite[datcite$`Document Type`=="Article" | datcite$`Document Type`=="Article; Early Access" | datcite$`Document Type`=="Letter" ,1:67]
  fulldat<-rbind(fulldat,dat1)
  dat1<-datcite[datcite$`Document Type`=="Article; Book Chapter" ,1:67]
  books<-rbind(books,dat1)
}
cat("Total of ",dim(fulldat)[1],"citations\n")
books<-books[!duplicated(books$`UT (Unique WOS ID)`),]
books$`Source Title`

# Removing duplicates
fulldat<-fulldat[!duplicated(fulldat$`UT (Unique WOS ID)`),]
  id1<-grep("Comets",fulldat$Authors)
  id2<-grep("Mentr",fulldat$Authors)
  id<-unique(c(id1,id2))
cat("Number of unique articles citing any of the",length(nampub),"articles:",dim(fulldat)[1],"including",length(id),"auto-citations \n")

removecol<-c(3:5,7,11:12,15:19,30,48:51,56:57)

write.table(fulldat[,-c(removecol)], file.path("uniquepublications.csv"), quote=TRUE, row.names=F)

# Annotated file

workdat<-read.table(file.path("uniquepublications_annotated.csv"),stringsAsFactors = FALSE, header=TRUE)

# Removing our 5 articles
cite.pubmedid<-c(18215437,0,20033477,16906454,27884052) #pubmedID for the 5 articles, one (Escolano) missing because first in series
cite.pubmedid<-cite.pubmedid[cite.pubmedid>0]
workdat1<-workdat[!(workdat$Pubmed.Id %in% cite.pubmedid),] #510 publications

```

### Extracting addresses of authors

```{r authorAddress}
# Extracting country of corresponding author
address<-workdat1$Reprint.Addresses
address<-strsplit(address,",",fixed=TRUE)
corcountry<-c()
for(i in address) {
  x<-i[length(i)]
  x<-unlist(strsplit(x,".",fixed=TRUE))
  corcountry<-c(corcountry,x[x!=""])
}
# Removing postal code from the US addresses
idx<-grep("USA",corcountry)
region<-corcountry
for(i in idx) {
  x<-unlist(strsplit(corcountry[i]," ",fixed=TRUE))
  x<-x[x!=""]
  if(length(x)==3 & !is.na(as.double(x[2]))) x<-x[-c(2)]
  if(length(x)!=2) cat("Problem with",corcountry[i])
  corcountry[i]<-paste(x[1],x[2])
  region[i]<-state.name[state.abb==x[1]]
}

table(corcountry)
corcountry2<-corcountry
corcountry2[grep("USA",corcountry)]<-"USA"
x<-table(corcountry2)
addr<-as.data.frame(x)
colnames(addr)<-c("country","nb")
```

## Plots for paper

### World map of the papers citing the foundation papers

Creating a world map and colouring the countries according to the number of publications

```{r worldMap}
# World map
worldmap<-map_data("world")
worldmap2<-worldmap[worldmap$lat>(-51),]

trim.leading <- function (x)  sub("^\\s+", "", x)
addr$country<-trim.leading(addr$country)

# Checking we have all regions
pmatch(addr$country,worldmap$region)
addr$country[addr$country=="Bosnia & Herceg"]<-"Bosnia and Herzegovina"
addr$country[addr$country=="Peoples R China"]<-"China"
addr$country[addr$country=="England"]<-"UK"
# Regroup Scotland and the UK
addr$nb[addr$country=="UK"]<-addr$nb[addr$country=="UK"]+addr$nb[addr$country=="Scotland"]
addr<-addr[addr$country!="Scotland",]
pmatch(addr$country,worldmap$region)

# Adding colours to the dataset... Continuous scale
npaper<-addr$nb[match(worldmap2$region,addr$country)]
npaper[is.na(npaper)]<-0
worldmap3<-cbind(worldmap2,npaper=npaper)

#ggplot( data=worldmap3, aes(x=long, y=lat, fill=npaper, group=group)) + geom_polygon(color="black", size=0.5) + scale_fill_gradient()

# Adding colours to the paper, regrouping in categories
npaper[npaper<5 & npaper>1]<-1
npaper[npaper<11 & npaper>=5]<-2
npaper[npaper<50 & npaper>=11]<-3
npaper[npaper>=50]<-4
npaper<-npaper+1
worldmap3<-cbind(worldmap3,ntab=npaper)

#cols<-c(colorRampPalette(c("white", "royalblue"))(6))
#cols<-cols[-c(1)]
mycolours<-brewer.pal(5,"Blues")

mymap<-ggplot( data=worldmap3, aes(x=long, y=lat, fill=as.factor(ntab), group=group)) + geom_polygon(color="black", size=0.5) + xlab("Longitude") + ylab("Latitude") + theme(legend.position = "top")  + scale_fill_manual(values = mycolours, name = "Number of articles",labels=c("None","1-4","5-10","11-49","50 or more"))

print(mymap)
```

### Histogram of number of publications per year

Histograms of number of (unique) publications citing at least one of the foundation articles (including autocitations, using workdat1 below will exclude the 5 autocitations):

```{r histPubnb}
histTime<-ggplot(data=workdat, aes(x=Publication.Year, fill=as.factor(Methodology.Application))) + geom_bar() + xlab("Publication Year") + ylab("Number of articles") +
    theme(legend.position = "top")  + scale_fill_manual(values=mycolours[c(2,5)],name = "Type of article",labels=c("Application","Methodology"))

print(histTime)

cat("Articles citing npde including the foundation articles, methodology versus applications:\n")
print(table(workdat$Methodology.Application))

cat("Articles citing npde (excluding the foundation articles), methodology versus applications:\n")
print(table(workdat1$Methodology.Application))
```

### Word cloud of keywords

Only on the applied articles. Word cloud to view :

- types of models PK or PD
- types of studies
- age groups: children/elderly
- ...

#### Harmonising keywords

Regrouping keywords in broad(er) categories, harmonising.

```{r keywordsData}
# Focus on applications using npde
workdat2<-workdat[workdat$Methodology.Application=="A",]
keyw<-c()
for(i in 1:dim(workdat2)[1]) {
  vec<-c(workdat2$Keywords.Plus[i],workdat2$Author.Keywords[i])
  vec<-vec[!is.na(vec)]
  if(length(vec)>1) vec<-paste(vec,collapse=";")
  vec<-unlist(strsplit(vec,";",fixed=TRUE))
  vec<-trimws(tolower(vec))
  keyw<-c(keyw,unique(vec))
}
# sort(unique(keyw[grep("pharmacodynamic",keyw)]))

# Table of keywords
keyw[grep("pk",keyw)]
keyw<-gsub("pharmacokinetics", "pk", keyw)
keyw<-gsub("pharmacokinetic", "pk", keyw)
keyw<-gsub("pharmacodynamics", "pd", keyw)
keyw<-gsub("pharmacodynamic", "pd", keyw)
keyw<-gsub("pkpd", "pk/pd", keyw)
keyw<-gsub("pk pd", "pk/pd", keyw)
keyw<-gsub("pk-pd", "pk/pd", keyw)
keyw<-gsub("modeling", "modelling", keyw)
keyw<-gsub("pharmacokinetick-pharmacodynamic", "pk/pd", keyw)
keyw<-gsub("pk/pd, pharmacokinetic-pharmacodynamic", "pk/pd", keyw)
keyw<-gsub("modeling", "model", keyw)
keyw<-gsub("modelling", "model", keyw)
keyw<-gsub("models", "model", keyw)
keyw<-gsub("healthy-volunteers", "healthy volunteers", keyw)
# Regrouping children, neonates, infants and preterm
keyw<-gsub("children", "paediatrics", keyw)
keyw<-gsub("neonates", "paediatrics", keyw)
keyw<-gsub("neonate", "paediatrics", keyw)
keyw<-gsub("pediatric-patients", "paediatrics", keyw)
keyw<-gsub("pediatrics", "paediatrics", keyw)
keyw<-gsub("pediatric", "paediatrics", keyw)
keyw<-gsub("preterm", "paediatrics", keyw)
keyw<-gsub("infants", "paediatrics", keyw)
keyw<-gsub("paediatrics", "paediatric", keyw)
keyw<-gsub("model based", "model-based", keyw)
keyw<-gsub("hiv-1", "hiv", keyw)
keyw<-gsub("infections", "infection", keyw)

# adults/children
keyw[keyw=="premenopausal women" | keyw== "african women" | keyw=="postmenopausal women"]<-"women"
keyw<-gsub("young-adults", "adult", keyw)
keyw<-gsub("adults", "adult", keyw)
keyw<-gsub("adult patients", "adult", keyw)
keyw<-gsub("adult chinese patients", "adult", keyw)
keyw<-gsub("chinese adult patients", "adult", keyw)
keyw[keyw== "paediatric pk population model"]<- "paediatric population pk model"
keyw[keyw=="chinese adult"]<-"adult"
keyw[keyw=="young paediatric" | keyw=="young-paediatric" | keyw=="zambian paediatric" | keyw=="paediatric paediatric" |keyw=="thai paediatric" | keyw=="paediatric inpatients" | keyw=="paediatric population" | keyw=="premature-paediatric" | keyw== "chinese premature paediatric" |keyw== "paediatric and young paediatric" | keyw=="younger patients" ]<-"paediatric"
keyw[keyw=="hiv-infected paediatric" | keyw=="hiv-infected children"]<-"hiv-infected paediatric"  
keyw[keyw=="paediatric and young paediatric" | keyw=="paediatric younger" | keyw=="premature paediatric"  |keyw=="chinese paediatric"   ]<-"paediatric"
keyw[keyw=="critically ill paediatric"]<-"critically-ill paediatric"
keyw[keyw=="newborn-paediatric" | keyw=="paediatric and adolescents" | keyw=="uninfected paediatric"]<-"paediatric"
keyw[keyw=="elderly hospitalized-patients" | keyw=="elderly-patients" | keyw=="older patients" | keyw=="older men"]<-"elderly"
keyw[keyw=="volunteers" | keyw=="healthy-subjects"]<-"healthy volunteers"

# ICU
keyw<-gsub("intensive care unit", "intensive care", keyw)
keyw<-gsub("intensive-care-unit", "intensive care", keyw)
keyw<-gsub("intensive-care patients", "intensive care", keyw)
keyw<-gsub("intensive care unit patients", "intensive care", keyw)
keyw<-gsub("intensive care patients", "intensive care", keyw)
keyw<-gsub("intensive-care", "intensive care", keyw)
keyw<-gsub("critically-ill patients", "intensive care", keyw)

keyw[keyw=="drug monitoring"]<-"therapeutic drug monitoring"
keyw[keyw=="drug-monitoring service"]<-"therapeutic drug monitoring"
keyw[keyw=="drug metabolism"]<-"drug-metabolism"
keyw[keyw=="drug development"]<-"drug-development"
keyw[keyw=="drug disposition"]<-"drug-disposition"
keyw[keyw=="drug interaction"| keyw=="drug interactions" | keyw=="drug-interactions" | keyw=="drug-drug interactions" | keyw=="drug-drug interaction"]<-"drug-interaction"
keyw[keyw=="drug"]<-"drugs"
keyw[keyw=="mediated drug disposition" | keyw=="target-mediated drug disposition" | keyw=="tmdd, target-mediated drug disposition"]<-"TMDD"

# PK/PD/model
keyw[keyw=="pk and pd"]<-"pk/pd"
keyw[keyw=="pd pd"]<-"pd"
keyw[keyw=="pk/pd (pk/pd) model"]<-"pk/pd model"
keyw[keyw=="nlme model"]<-"nonlinear mixed-effect model"
keyw[keyw=="pharmacodynamic model" | keyw=="pd model"| keyw=="pharmacodynamic model"]<-"pd model"
keyw[keyw=="mixed effects model" | keyw=="mixed-effects model" | keyw=="mixed effect model" | keyw=="mixed-effect model"]<-"mixed-effects model"
keyw[keyw=="nonlinear mixed effects model" | keyw=="nonlinear mixed-effect model" | keyw=="nonlinear mixed effect model" | keyw=="nonlinear mixed-effects model" | keyw=="nonlinear mixed effects model (nonmem)" | keyw=="nonlinear mixed model" | keyw=="nonlinear mixed-effects model"]<-"nonlinear mixed-effects model"
keyw[keyw=="pd model"]<-"pd model"
keyw[keyw=="pk/pd model"]<-"pk/pd model"
keyw[keyw=="pk/pd, pk/pd"]<-"pk/pd"
keyw[keyw=="population pk and pd"]<-"population pk/pd"
keyw[keyw=="kpd model" | keyw=="k-pd model"]<-"k-pd model"
keyw[keyw=="pk/pd analysis"]<-"pk/pd"
keyw[keyw=="population pk analysis"]<-"population pk"
keyw[keyw=="population pk model analysis" ]<-"population pk model"
keyw[keyw=="pk/pd relationship"]<-"pk/pd model"
keyw[keyw=="probability of pk/pd target attainment" | keyw=="pk/pd target attainment rate"]<-"pk/pd target attainment"

# Transplantation
keyw[keyw=="renal transplantation" | keyw=="renal-transplantation" | keyw=="renal-transplant recipients" | keyw=="kidney transplantation" ]<-"renal transplant"
keyw[keyw=="stem cell transplantations"]<-"stem-cell transplantation" 
keyw[keyw=="liver-transplantation" | keyw=="liver transplant"]<-"liver transplantation"
keyw[keyw=="marrow-transplantation"]<-"marrow transplantation"
keyw[keyw=="paediatric liver-transplantation" | keyw=="paediatric liver-transplantation" | keyw=="paediatric liver transplantation"]<-"paediatric liver transplantation"
keyw[keyw=="marrow-transplantation"]<-"marrow transplantation"
keyw[keyw=="kidney-transplant recipients" | keyw=="kidney-transplantation" | keyw=="kidney transplant patients"]<- "kidney transplantation"
keyw[keyw=="hematopoeitic stem cell transplantation" | keyw=="hematopoietic-cell transplantation" | keyw=="bone-marrow-transplantation"| keyw=="marrow transplantation" |keyw=="stem-cell transplantation"]<-"hematopoietic stem cell transplantation"
keyw[keyw=="transplant recipients" | keyw=="transplant patients"]<-"transplantation"

# Cancer
keyw[keyw=="breast-cancer patients" | keyw=="breast-cancer"| keyw== "her2-positive breast-cancer" | keyw=="previously untreated metastatic breast cancer" | keyw=="metastatic breast-cancer" ]<-"breast cancer"
keyw[keyw=="cancer patients" | keyw=="cancer-patients" | keyw=="cancer clinical-trials"]<-"cancer"
keyw[keyw=="lung-cancer patients" | keyw=="lung-cancer" | keyw=="cell lung-cancer"]<-"lung cancer"
keyw[keyw=="colorectal-cancer patients" | keyw== "colorectal-cancer"|keyw=="metastatic colorectal-cancer" ]<-"colorectal cancer"
keyw[keyw=="anti-cancer agents" | keyw=="anticancer agents" ]<-"anticancer drugs"
keyw[keyw=="ovarian epithelial cancer" | keyw=="ovarian-cancer patients" |keyw=="ovarian-cancer"]<-"ovarian cancer" 

# Animals
keyw[keyw=="murine model" | keyw=="juvenile mice"]<-"mice"
keyw[keyw=="animal-model" | keyw=="animals"]<-"animal"

keyw[keyw=="medication" | keyw=="medications"]<-"drugs"

# Parameters
keyw[keyw=="drug clearance" | keyw=="plasma-clearance" | keyw=="renal drug clearance" | keyw=="renal clearance" | keyw=="clearance rate" | keyw=="oral clearance" | keyw=="apparent clearance" | keyw=="plasma-clearance"]<-"clearance"
keyw[keyw=="drug-metabolism"| keyw=="first-pass metabolism"]<-"metabolism"                 
keyw[keyw=="drug-disposition" | keyw=="disposition"]<-"elimination"
keyw[keyw=="drug absorption"|keyw=="oral-drug absorption"|keyw=="systemic absorption"]<- "absorption"           
keyw[keyw=="glomerular filtration rate (gfr)" | keyw=="glomerular-filtration-rate" | keyw=="glomerular rate filtering" | keyw=="glomerular-filtration" | keyw=="glomerular filtration"]<-"glomerular filtration rate"
keyw[keyw=="severe sepsis" | keyw=="late-onset sepsis"]<-"sepsis"
keyw[keyw=="human-plasma" | keyw=="human plasma"]<-"plasma"

# Model purpose
keyw[keyw=="model validation" ]<-"model evaluation" 
keyw[keyw=="mathematical-model"]<-"mathematical model"
keyw[keyw=="mechanism-based model" ]<-"mechanistic model" 
keyw[keyw=="single-dose pk"]<-"single dose pk"
keyw[keyw=="vitro" | keyw=="in vitro experiments" | keyw=="in-vitro activity" | keyw=="in-vitro" | keyw=="in-vitro metabolism" | keyw=="vitro dynamic-model"]<-"in vitro"
keyw[keyw=="pharmacogenetic" | keyw=="pk-pharmacogenetics"]<-"pharmacogenetics"
keyw[keyw=="genetic-polymorphism"]<-"genetic polymorphisms"

# Study
keyw[keyw=="drug dose-response relationship" | keyw=="dose-response relationship" ] <- "dose-response"
keyw[keyw=="clinical-pharmacology"]<-"clinical pharmacology"
keyw[keyw=="clinical trials" | keyw=="clinical-trials"]<-"clinical trial"
keyw[keyw=="phase i clinical trial" | keyw=="phase-i" | keyw=="phase-ii" | keyw=="phase i" | keyw=="phase-i trial" | keyw=="phase-iii trial"]<-"clinical trial"
keyw[keyw=="randomized-controlled-trial" | keyw=="randomized controlled-trial" | keyw=="randomized-trial" | keyw=="placebo-controlled trial" | keyw=="cancer clinical-trials" |
keyw=="controlled-trial"]<-"randomized controlled trial"
keyw[keyw=="trials"]<-"trial"
keyw[keyw=="clinical-efficacy"]<-"efficacy"
keyw[keyw=="clinical-outcomes" | keyw=="clinical-response"]<-"clinical response"
keyw[keyw=="drug-development" | keyw=="model-based drug development"]<-"drug development"
# Dose individualisation
keyw[keyw %in% c("dose individualisation","dose adjustment","adaptive dosing","dosage optimization", "dose individualization","dosage adaptation","dosage prediction","individualized dosing","optimal dose","optimal dosing","dosage adjustment","dose optimization","individualized drug dosing","dosing optimization")]<-"dose individualisation"
# Dose selection (population)
keyw[keyw %in% c("drug dosage calculations","dosage","dosing recommendation","child appropriate dosage forms","dosing chart","dosage recommendations","dose rationale","dose requirements","dosing","dosing guideline","dosing guidelines","dosing rational","dosing regimen","dosing regimens","dosing schedule","dosing scheme","dosing strategy","dosing table" ,"drug dosage calculations","drug dosing","initial dosage regimens","recommended dose")]<-"drug selection"

keyw[keyw=="safety pharmacology"|keyw=="long term safety"]<-"safety"
keyw[grep("toxicity",keyw)]<-"toxicity"

# Conditions
keyw[keyw=="migraine disorders" | keyw=="migraine prevention"]<-"migraine"
keyw[keyw=="hiv-infected adults patients" | keyw=="hiv-infected individuals" | keyw=="immunodeficiency-virus-infection"]<-"hiv-infected adult"
keyw[keyw=="infectious diseases" | keyw=="infected patients" ]<-"infection"
keyw[keyw=="virus-infected patients" | keyw=="virus-infected subjects" | keyw=="virus-infection"  ]<-"viral infection" # technically, hiv, influenza, hepatitis could also go in viral infection
keyw[keyw=="impaired renal-function" | keyw=="renal-function"]<-"renal function"
keyw[keyw=="renal dysfunction" | keyw== "impaired renal-function" | keyw=="chronic-renal-failure" | keyw=="renal-failure"]<-"renal impairment"
keyw[keyw=="human liver"]<-"liver"
keyw[keyw=="chronic hepatitis-b" | keyw=="chronic hepatitis-c" | keyw=="hepatitis c infection" | keyw=="hepatitis c virus (hcv)"]<-"hepatitis"
keyw[keyw=="morbid obesity" | keyw=="obese-patients" | keyw=="morbidly obese-patients"]<-"obesity"
keyw[keyw=="overweight"]<-"under/overweight"
keyw[keyw=="crohn's disease" | keyw=="crohns-disease"]<-"crohn disease"
keyw[keyw=="inflammatory-bowel-disease"]<-"inflammatory bowel disease"
keyw[keyw=="parkinson's disease" | keyw=="parkinsons-disease"]<-"Parkinson's disease"

# Treatments
keyw[keyw=="clinical-evaluation"]<-"clinical evaluation"
keyw[keyw=="multidrug-resistant" | keyw=="multidrug-resistance" | keyw=="mrp2 (multi-drug resistance protein-2)/abcc2" | keyw=="multidrug-resistance-associated proteins" | keyw== "drug resistance"]<-"drug resistance"
keyw[keyw=="antituberculosis drugs"]<-"antitubercular drugs"
keyw[keyw=="antibody-drug conjugate"]<-"antibody drug conjugate"                     
keyw[keyw=="plasma clozapine"]<-"clozapine"
keyw[keyw=="lopinavir plasma-concentrations"]<-"lopinavir"
keyw[keyw=="combinations" | keyw=="combination therapy" | keyw=="fixed-dose combination"]<-"combination"
keyw[keyw=="recombinant-human-erythropoietin"]<-"recombinant erythropoietin"
keyw[keyw=="antimicrobial agents" | keyw=="combination antimicrobial therapy" | keyw=="nebulized antimicrobial agents"]<-"antimicrobial therapy"
keyw[keyw=="microbial response"]<-"antimicrobial activity"

# Other
keyw[grep("polymorphism",keyw)]<-"genetic polymorphism"
keyw[grep("pain",keyw)]<-"pain"
keyw[keyw=="p-glycoprotein inhibition"]<-"p-glycoprotein inhibitor"
keyw[keyw %in% c("developmental changes","developmental expression","developmental pharmacology")]<-"developmental pk"
keyw[keyw %in% c("monte carlo simulation","monte-carlo simulation","monte-carlo-simulation","monte carlo simulations")]<-"monte carlo simulation"
keyw[keyw=="simulations" | keyw=="computer simulation" ]<-"simulation"
keyw[keyw %in% c("serum concentrations","serum concentration","plasma-concentrations", "blood-concentrations", "unbound serum concentrations","target site concentrations","free drug concentration","effect-site concentration","drug concentrations")]<-"concentrations"
keyw[keyw=="trough concentrations"]<-"trough concentration"
  
if(FALSE) {
keyw[keyw==""]<-""
}

keyw[keyw=="pk"]<-"pharmacokinetics"
keyw[keyw=="pd"]<-"pharmacodynamics"


# Removing duplicate keyword
keyw<-keyw[keyw!="1-infected paediatric"]

```

```{r keywordsCheck}
# Looking at some specific keywords
sort(unique(keyw[grep("pd",keyw)]))
sort(unique(keyw[grep("pk",keyw)]))
#sort(unique(keyw[grep("drug monitoring",keyw)]))
sort(unique(keyw[grep("adult",keyw)])) # regrouped and split
sort(unique(keyw[grep("paediatric",keyw)])) # regrouped and split
sort(unique(keyw[grep("neonate",keyw)]))
sort(unique(keyw[grep("transplant",keyw)]))
sort(unique(keyw[grep("cancer",keyw)]))
sort(unique(keyw[grep("tumo",keyw)]))
sort(unique(keyw[grep("development",keyw)]))
sort(unique(keyw[grep("dog",keyw)])) # no cat
sort(unique(keyw[grep("mice",keyw)])) # no mouse
sort(unique(keyw[grep("monkey",keyw)])) # no primate, monkey
sort(unique(keyw[grep("monitoring",keyw)]))

sort(unique(keyw[grep("evaluation",keyw)]))
sort(unique(keyw[grep("internal",keyw)]))
sort(unique(keyw[grep("intensive",keyw)]))
sort(unique(keyw[grep("drug",keyw)]))
sort(unique(keyw[grep("model",keyw)]))
sort(unique(keyw[grep("distribution",keyw)]))
sort(unique(keyw[grep("volunteer",keyw)]))
sort(unique(keyw[grep("recipient",keyw)]))

```

#### Keywords covering several notions

Some keywords cover different notions. The usual text parsing approach would be to cut strings into single words, but here trying to keep expressions => split the keywords into 2

```{r keywordsSplit}
# Splitting keywords
splitkw<-function(vec, pattern, replist) {
  idx<-which(vec==pattern)
  vec[idx]<-replist[1]
  for(i in 2:length(replist))
    vec<-c(vec,rep(replist[i],length(idx)))
  return(vec)
}
# splitkw(1:5,2,c(10,11,12)) # Testing the function
# splitkw(rep(1:5, each=3),2,c(10,11,12))

# 
keyw<-gsub("models", "model", keyw)
keyw<-splitkw(keyw,"population pk/pd model",c("population pk/pd","model"))
keyw<-splitkw(keyw,"population pk model",c("population pk","model"))
keyw<-splitkw(keyw,"pk/pd models",c("pk/pd","model"))
keyw<-splitkw(keyw,"pk models",c("pk","model"))
keyw<-splitkw(keyw,"pk model",c("pk","model"))
keyw<-splitkw(keyw,"pbpk model",c("pbpk","model"))
keyw<-splitkw(keyw,"pd models",c("pd","model"))
keyw<-splitkw(keyw,"whole-body physiologically based pk",c("pbpk","model"))

# Adult + specific disease
keyw<-splitkw(keyw,"hiv-infected adult",c("hiv","adult","infectious disease"))
keyw<-splitkw(keyw,"adult t-cell lymphoma",c("adult","cancer"))
keyw<-splitkw(keyw,"infected women",c("women","infection"))
# Children + specific condition
keyw<-splitkw(keyw,"paediatric intensive care",c("paediatric","intensive care"))
keyw<-splitkw(keyw,"paediatric population pk model",c("paediatric","population pk","model"))
keyw<-splitkw(keyw,"hiv-infected paediatric",c("hiv","paediatric"))
keyw<-splitkw(keyw,"paediatric intensive care",c("paediatric","intensive care"))
keyw<-splitkw(keyw,"critically ill paediatric",c("paediatric","intensive care"))
keyw<-splitkw(keyw,"critically-ill paediatric",c("paediatric","intensive care"))
keyw<-splitkw(keyw,"paediatric liver transplantation",c("paediatric","liver transplant","transplantation"))
keyw<-splitkw(keyw,"paediatric oncology patients" ,c("paediatric","cancer"))
keyw<-splitkw(keyw, "paediatric brain tumors",c( "paediatric","cancer"))
keyw<-splitkw(keyw,"paediatric pharmacology",c("paediatric","pharmacology"))
keyw<-splitkw(keyw,"paediatric clinical pharmacology",c("paediatric","pharmacology"))
keyw<-splitkw(keyw,"birth-weight paediatric",c("paediatric","neonates"))
keyw<-splitkw(keyw,"breast-fed paediatric",c("paediatric","neonates"))
keyw<-splitkw(keyw,"paediatric pharmacology",c("paediatric","pharmacology"))
keyw<-splitkw(keyw,"paediatric clinical pharmacology",c("paediatric","pharmacology"))
keyw<-splitkw(keyw,"paediatric anesthesia",c("paediatric","anesthesia"))
keyw<-splitkw(keyw,"paediatric drug development",c("paediatric","drug development"))
keyw<-splitkw(keyw,"paediatric population pk model",c("paediatric","population pk","model"))
keyw<-splitkw(keyw,"paediatric liver transplantation",c("paediatric","liver transplantation"))
keyw<-splitkw(keyw,"paediatric migraineurs",c("paediatric","migraine"))
keyw<-splitkw(keyw,"childhood migraine",c("paediatric","migraine"))
keyw<-splitkw(keyw,"neonatal sepsis",c("paediatric","neonates","sepsis"))

# Animals
keyw<-splitkw(keyw,"mice",c("mice","animal"))
keyw<-splitkw(keyw,"veterinary drugs",c("animal","drugs"))

# Cancer

# Transplantation (any)
keyw<-splitkw(keyw,"renal transplant",c("renal transplant","transplantation"))
keyw<-splitkw(keyw,"adult renal transplant recipients",c("renal transplant","transplantation","adult"))
keyw<-splitkw(keyw,"heart-transplantation",c("heart-transplantation","transplantation"))
keyw<-splitkw(keyw,"liver transplantation",c("liver transplantation","transplantation"))
keyw<-splitkw(keyw,"hematopoietic stem cell transplantation",c("hematopoietic stem cell transplantation","transplantation"))
keyw<-splitkw(keyw,"solid-organ transplantation",c("solid-organ transplantation","transplantation"))

# Drugs
keyw<-splitkw(keyw,"piperacillin population pk",c("piperacillin","population pk"))
keyw<-splitkw(keyw,"gentamicin pk",c("gentamicin","pk"))
keyw<-splitkw(keyw,"anticancer drug resistance",c("cancer","drug resistance"))

#keyw<-splitkw(keyw,"",c())

if(FALSE) {
keyw<-splitkw(keyw,"",c())
keyw<-splitkw(keyw,"",c())
keyw<-splitkw(keyw,"",c())
}
```

#### Creating word cloud figure

```{r keywordsCloud.tweaked}
dm<-data.frame(keyw=unique(keyw),count=NA,stringsAsFactors = FALSE)
for(i in unique(keyw)) {
  dm$count[dm$keyw==i]<-length(keyw[keyw==i])
}
dm<-dm[order(dm$count, decreasing = TRUE),]
d <- data.frame(word = dm$keyw, freq=dm$count)
d[d$freq>=10,]
cat(dim(d[d$freq>=10,])[1],"keywords occurring 10 times or more, and", dim(d[d$freq<10 & d$freq>=5,])[1],"keywords occurring between 5 and 10 times\n")

#d[d$freq<=10 & d$freq>5,]

set.seed(1234567)

wordcloud(words = d$word, freq = d$freq, min.freq = 5,scale=c(2,.5), 
    max.words=50, random.order=FALSE, rot.per=0.35, colors=brewer.pal(8, "Dark2"))
```


## Diagnostic plots for viral load example

### Install package in development mode

This is a snapshot of the version of npde released at the time of the revision of the paper submitted to AAPS Journal. The npde package version 3.0 has been submitted to CRAN and should be available shortly on CRAN. We use development mode to avoid messing up any existing installation of npde.

```{r devModeInstall}
dev_mode() # development mode
install.packages(pkgs="npde_3.0.tar.gz",repos=NULL) # local installation
# install.package("npde) # installation from CRAN
library(npde)
```

### Default diagnostic plots for npd for the viral load example used in the paper

```{r virload, echo=FALSE}
#Viral load data
data(virload)
data(simvirload)

#Plotting the data
par(mfrow=c(1,2))
plot(Log_VL~Time,data=virload,xlab="Time (d)",ylab="Viral loads, base 10 log-scale (cp/mL)")

# Calling autonpde with dataframes
yvir<-autonpde(virload,simvirload,iid=1,ix=2,iy=3,boolsave=FALSE, units=list("hr","cp/mL, log 10 base"))

plot(yvir)

# Saving figures to folder figsDir
figsDir<-getwd() #figsDir<-"/home/eco/xtex/npde/aaps20/figs"

if(FALSE) {
  namfile<-file.path(figsDir,paste("npdeDefault.","eps",sep=""))
  cairo_ps(file = namfile, width = 11, height = 8, pointsize=20)
  plot(yvir)
  dev.off()
}
```

### npd for the warfarin example in the documentation

Running the example:

```{r warfarin, echo=FALSE}
# Warfarin with covariate model
data(warfarin)
data(simwarfarinCov)
wcov<-autonpde(namobs=warfarin,namsim=simwarfarinCov, iid=1,ix=2,iy=4,icov=c(3,6:8),namsav="warfCov", units=list(x="hr",y="mg/L", covariates=c("mg","kg","-","yr")))
```

Reference plot for subject 2

```{r warfarinRefPlot, echo=FALSE}
plot.tnpde<-plot(wcov, plot.type="x.scatter", ref.prof=list(id=2), main="tnpd with reference profile ID=2")
plot.vpc<-plot(wcov, plot.type="vpc", main="VPC")
grid.arrange(grobs=list(plot.tnpde, plot.vpc), nrow=1, ncol=2)

if(FALSE) {
  namfile<-file.path(figsDir,paste("referenceProfile.","eps",sep=""))
  cairo_ps(file = namfile, width = 11, height = 8, pointsize=20)
  grid.arrange(grobs=list(plot.tnpde, plot.vpc), nrow=1, ncol=2)
  dev.off()
}
```

The code below shows several diagnostic plots for covariates:

- plots like scatterplots or distribution plots can be split over categories of covariates
- we can also plot the distribution of eg npd versus the categories as boxplots

```{r warfarinCovPlot, echo=FALSE}
xwt.scatt<-plot(wcov, plot.type="x.scatter", covsplit=TRUE, which.cov="wt")
xwt.qqplot<-plot(wcov, plot.type="qqplot", covsplit=TRUE, which.cov="wt")
xwt.box<-plot(wcov, plot.type="covariates", which.cov="wt")
xsex.scatt<-plot(wcov, plot.type="x.scatter", covsplit=TRUE, which.cov="sex")
xsex.qqplot<-plot(wcov, plot.type="qqplot", covsplit=TRUE, which.cov="sex")
xsex.box<-plot(wcov, plot.type="covariates", which.cov="sex")

# Sometimes object is returned as a list, if the message "only 'grobs' allowed in "gList"" appears then try using [[1]] after the object
grid.arrange(grobs=list(xwt.scatt, xsex.scatt, xwt.qqplot, xsex.qqplot, xwt.box[[1]],  xsex.box[[1]]), nrow=3, ncol=2)

if(FALSE) {
  namfile<-file.path(figsDir,paste("covariateplots.","eps",sep=""))
  cairo_ps(file = namfile, width = 11, height = 8, pointsize=20)
  grid.arrange(grobs=list(xwt.scatt, xsex.scatt, xwt.qqplot, xsex.qqplot, xwt.box[[1]],  xsex.box[[1]]), nrow=3, ncol=2)
  dev.off()
}

# grid.arrange(grobs=list(xwt.scatt, xwt.qqplot, xwt.box[[1]], xsex.scatt, xsex.qqplot, xsex.box[[1]]), nrow=2, ncol=3)
# grid.arrange(grobs=list(xwt.scatt, xwt.qqplot, xsex.scatt, xsex.qqplot), nrow=2, ncol=2)
```

## End of code
Exit development mode.

```{r devModeStop}
dev_mode() # development mode
```
