Nothing
#Copyright (c) 2009-2022 Sebastien Bihorel
#All rights reserved.
#
#This file is part of scaRabee.
#
# scaRabee is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# scaRabee is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with scaRabee. If not, see <http://www.gnu.org/licenses/>.
#
estimation.plot <- function(problem=NULL,Fit=NULL,files=NULL){
# Set some default variable value
nestparam <- length(Fit$estimations)
# Create ... vs iteration plots
# Determine the most esthetical arrangement of plots based on nestparam
tmplayout <- get.layout(nplot=nestparam + 1)
# Import and reshape data to be plotted
data <- read.table(file=files$iter,
header=TRUE,
as.is=TRUE,
sep=',',
check.names=FALSE)
est.names <- problem$init$names[which(problem$init$isfix==0)]
data <- reshape(data[,c('Objective function',est.names,'ID')],
idvar='iteration',
times=c('Objective function',est.names),
timevar='variable',
varying=list(c('Objective function',est.names)),
direction='long')
data$variable <- factor(data$variable,levels=unique(data$variable))
firstiterations <- data$iteration[match(unique(data$ID),data$ID)]
for (id in unique(data$ID)){
data$iteration[which(data$ID==id)] <- data$iteration[which(data$ID==id)] -
firstiterations[id]
}
names(data)[3] <- 'value'
# Objective function and parameters vs iteration plot
pdf(file='01.Iterations.pdf',
width=6.5,
height=9)
for (i in unique(data$ID)) {
if (problem$method=='subject'){
ititle <- paste('Subject', i)
} else {
ititle <- 'Population'
}
itplot <- xyplot(value~iteration|variable,
data=data,
subset=which(data$ID==i),
as.table=TRUE,
type='l',
col=4,
scales=list(x=list(relation='free'),
y=list(relation='free')),
strip=strip.custom(style=1,bg=0),
xlab='Iteration',
ylab='Value / Estimate',
main=paste('Gradient monitoring - ',ititle,sep=''),
layout=tmplayout)
print(itplot)
}
dev.off()
# Create diagnostic plots
# Import data to be plotted
data <- read.table(file=files$pred,
header=TRUE,
as.is=TRUE,
sep=',',
check.names=FALSE)
# Creates prediction plots
pdf(file='02.Predictions.pdf',
width=6.5,
height=9)
for (iid in unique(data$ID)){
subdata <- data[which(data$ID==iid),]
for (itrt in unique(subdata$TRT)){
if (problem$method=='subject'){
ititle <- paste('Subject', iid, '- Treatment', itrt)
} else {
ititle <- paste('Population - Treatment', itrt)
}
trtdata <- subdata[which(subdata$TRT==itrt),]
trtdata$DVID <- factor(trtdata$DVID,levels=unique(trtdata$DVID))
noutput <- nlevels(trtdata$DVID)
tmplayout <- get.layout(nplot=noutput)
myplot <- xyplot(DV+IPRED~TIME|DVID,
data=trtdata,
as.table=TRUE,
type=c('p','l'),
pch=c(3,NULL),
col=c(1,4),
distribute.type=TRUE,
scales=list(x=list(relation='free'),
y=list(relation='free')),
strip=strip.custom(var.name='Output',
style=1,
bg=0,
sep=' ',
strip.name=c(TRUE,TRUE)),
xlab='Time',
ylab='Observations / Predictions',
main=ititle,
layout=tmplayout)
print(myplot)
}
}
dev.off()
# Creates residuals plots
pdf(file='03.Residuals.pdf',
width=6.5,
height=9)
for (iid in unique(data$ID)){
subdata <- data[which(data$ID==iid),]
for (itrt in unique(subdata$TRT)){
trtdata <- subdata[which(subdata$TRT==itrt),]
for (icmt in unique(trtdata$DVID)){
if (problem$method=='subject'){
ititle <- paste('Subject', iid, '- Treatment', itrt, '- Output', icmt)
} else {
ititle <- paste('Population - Treatment', itrt, '- Output', icmt)
}
cmtdata <- trtdata[which(trtdata$DVID==icmt),]
rplot1 <- xyplot(IPRED~DV,
data=cmtdata,
panel=function(x,y,...){
panel.xyplot(x,y,type='p',pch=3,col=4)
panel.abline(c(0,1),type='l',col=2)
},
strip=NULL,
xlab='Observations',
ylab='Predictions')
rplot2 <- xyplot(WRES~TIME,
data=cmtdata,
panel=function(x,y,...){
panel.xyplot(x,y,type='p',pch=3,col=4)
panel.abline(c(0,0),type='l',col=2)
},
strip=NULL,
xlab='Time',
ylab='Weighted residuals')
rplot3 <- xyplot(WRES~DV,
data=cmtdata,
panel=function(x,y,...){
panel.xyplot(x,y,type='p',pch=3,col=4)
panel.abline(c(0,0),type='l',col=2)
},
strip=NULL,
xlab='Observations',
ylab='Weighted residuals')
rplot4 <- xyplot(WRES~IPRED,
data=cmtdata,
panel=function(x,y,...){
panel.xyplot(x,y,type='p',pch=3,col=4)
panel.abline(c(0,0),type='l',col=2)
},
strip=NULL,
xlab='Predictions',
ylab='Weighted residuals')
print(rplot1,position=c(0,0.1,0.5,0.5),more=TRUE)
grid.text(label=ititle,
x=0.5,y=0.95,just='center',gp=gpar(cex=1.5))
print(rplot2,position=c(0.5,0.1,1,0.5),more=TRUE)
print(rplot3,position=c(0,0.5,0.5,0.9),more=TRUE)
print(rplot4,position=c(0.5,0.5,1,0.9),more=FALSE)
}
}
}
dev.off()
# Print a message to screen
cat(sprintf('\n%s%s\n%s\n','Diagnostic figures have been created and saved in: ',
getwd(),'You may open and edit them at your convenience.'))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.