Este código de R corresponde a la figura 6.17zz del libro Análisis de datos con el programa estadístico R: una introducción aplicada de Salas-Eljatib (2021).

Esta figura corresponde a un panel de 9 gráficos de curvas con las siguientes características: (i) se diferencia cada columna con letras en negrita en la esquina superior derecha del gráfico; (ii) hay tres curvas por grafico con distintos colores y tipos; (iii) existe un recuadro con un zoom de las curvas en la otra esquina superior izquierda del gráfico; (iv) la etiqueta del eje-Y es a color; y (v) la etiqueta del eje-X utiliza símbolos estadísticos con letras griegas.

Cargando los datos de ejemplo

El estudio de Gregoire and Salas (2009) involucró investigar las propiedades estadisticas de tres estimadores de razón en el contexto de errores de medición. Cada columna del panel representa a una funcion de distribución de probabilidades especifica, mientras que cada fila del panel es una propiedad estadística de los estimadores. Además, en un recuadro se expresa una porción de las curvas pero ahora representando los estadísticos en términos porcentuales. Los datos necesarios para reproducir esta figura se encuentra en el objeto de datos measueratio.RData .

library(datana)
load("measueratio.RData")

El gráfico

Ahora se procede a realizar el gráfico.

#primero: un par de definiciones
lty.rat1=1;lty.rat2=6;lty.rat3=5
col.rat1="red";col.rat2="blue";col.rat3="black"
loc.line.lab.y=1.8 #2.2

label.x.rand=expression(paste(sigma [delta], " / ", sigma [x]))
mar.2=3.2;mar.3=2;mar.1=2;mar.4=0.4

#segundo: la produccion de la figura
oldpar <- par(mfrow=c(3,3), las=1,  oma = c(mar.1,0.2,mar.3,mar.4), mar = c(2, mar.2, 0.2, 0.2))


# Fig.8.(firs-row panel) - Bias ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
y.lim.min=1;y.lim.max=3.5#5
y.lim.min.inn=0;y.lim.max.inn=3
loc.line.lab.y=2.2 #2.4
pos.x.leg=0.29;pos.y.leg=y.lim.max*0.99;
pos.x.leg1=0.28;pos.y.leg1=y.lim.max*0.95;
lim.low.inn=8.7#12; the bigger the bigger the inner-plot
lim.left.inn=mar.2;lim.up.inn=0.25;lim.rig.inn=6
#abajo (> mas arriba), lado-izq,arriba,lado-der
per.pos.x.ylabel.inn=0.83;per.pos.y.ylabel.inn=0.93


# a. Uniform measurement error ~~~~~~~~~~~
plot(sddeltaUnifME/sigmax, BiasRatioUnifME.rat1, type = "l", xlim=c(0,0.3),                   
     xlab="", ylab="",  main="", col=col.rat1, ylim = c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaUnifME/sigmax, BiasRatioUnifME.rat2, type = "l",col = col.rat2, lty = lty.rat2)
#title(xlab=list(label.x.rand), col="black",line = 2, font=2, cex.lab=1.250)
title(ylab=list(expression("Raz\u00f3n de sesgo"), col="brown"), line =loc.line.lab.y, font=2, cex.lab=1.0)
text(pos.x.leg, pos.y.leg, expression(bold((a))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #0.83 for x



#b1=paste("Uniforme","                         Gaussian",sep="     ")
#c1=paste(b1       ,"                           Beta",sep="      ")
#mtext(c1,font=2,side=3,outer=TRUE)

# the outer top label -----------
b1=paste("                ", "Uniforme","                           Gauss",sep="       ")
c1=paste(b1       ,"                               Beta",sep="        ")

label.x.rand1=c1

mtext(label.x.rand1,font=2,side=3,outer=TRUE,adj=0)
# end of the ouer botton label -------


# Add an inner plot - Uniform ME
#opar <- par(new=T, oma=c(0, 0, 0, 0))
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn))                     
plot(sddeltaUnifME/sigmax, BiasPercentUnif_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))
#xaxs = "i", yaxs = "i"
#abline(h= c(0,1,2), col="black", lty = 3)
lines(sddeltaUnifME/sigmax, BiasPercentUnif_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
#title(ylab=list(expression("Bias (%)"),col="black"), line = 2, font=2, cex.lab=1.250)
par(opar)
# end of 'add an inner plot - Uniform ME'


# b. Normal measurement error ~~~~~~~~~~~
plot(sddeltaNormME/sigmax, BiasRatioNormME.rat1, type = "l", xlim=c(0,0.3),                   
     xlab="", ylab="",main="",col=col.rat1, ylim =c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaNormME/sigmax, BiasRatioNormME.rat2, type = "l", col = col.rat2, lty = lty.rat2)
text(pos.x.leg, pos.y.leg, expression(bold((b))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #0.83 for x

# Add an inner plot - Normal ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaNormME/sigmax, BiasPercentNorm_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)
lines(sddeltaNormME/sigmax, BiasPercentNorm_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)

par(opar)
# end of 'add an inner plot - Normal ME'



# c. Beta measurement error ~~~~~~~~~~~
plot(sddeltaBetaME/sigmax, BiasRatioBetaME.rat1, type = "l", xlim=c(0,0.3),                   
     xlab="", ylab="",  main="", col=col.rat1, ylim = c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaBetaME/sigmax, BiasRatioBetaME.rat2, type = "l", col = col.rat2, lty = lty.rat2)
text(pos.x.leg, pos.y.leg, expression(bold((c))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #0.83 for x


# Add an inner plot - Beta ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaBetaME/sigmax, BiasPercentBeta_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)
lines(sddeltaBetaME/sigmax, BiasPercentBeta_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)

par(opar)
# end of 'add an inner plot - Beta ME'
# end of 'Fig.8.1 - Bias' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# Fig.3.(second-row panel) - Std.Error ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
y.lim.min=1;y.lim.max=1.9#2
y.lim.min.inn=3;y.lim.max.inn=6
loc.line.lab.y=2.2 #2.4
pos.x.leg=0.29;pos.y.leg=y.lim.max*0.99;
pos.x.leg1=0.28;pos.y.leg1=y.lim.max*0.95;
lim.low.inn=8.7;#12
lim.left.inn=mar.2;lim.up.inn=0.2;lim.rig.inn=6
per.pos.x.ylabel.inn=0.83;per.pos.y.ylabel.inn=0.98



# a. Uniform measurement error ~~~~~~~~~~~
plot(sddeltaUnifME/sigmax, StdErrRatioUnifME.rat1,col=col.rat1,
     type = "l", lty=1, xlim=c(0,0.3), xlab="",ylab="",main="",ylim = c(y.lim.min,y.lim.max)) 
#abline(h=1, col="red", lty = 3) 

lines(sddeltaUnifME/sigmax, StdErrRatioUnifME.rat2, col=col.rat2,type = "l", lty=lty.rat2)    
lines(sddeltaUnifME/sigmax, StdErrRatioUnif_n7_ME_term1only.rat3,
      col=col.rat3,type = "l", lty=lty.rat3)

#title(xlab=list(label.x.rand), col="black",line = 2, font=2, cex.lab=1.250)
title(ylab=list(expression("Raz\u00f3n de Err.St."),col="brown"),
      line = loc.line.lab.y, font=2, cex.lab=1)
text(pos.x.leg, pos.y.leg, expression(bold((a))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #0.83 for x



# Add an inner plot - Uniform ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaUnifME/sigmax, StdErrPercentUnif_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1, ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))
lines(sddeltaUnifME/sigmax,StdErrPercentUnif_n7_ME.rat2, type = "l",
      lty =lty.rat2,xlab="",ylab="",col=col.rat2)
lines(sddeltaUnifME/sigmax,StdErrPercentUnif_n7_ME.rat3, type = "l",
      lty =lty.rat3,xlab="",ylab="",col=col.rat3)

par(opar)
# end of 'add an inner plot - Uniform ME'                             




# b. Normal measurement error ~~~~~~~~~~~
plot(sddeltaNormME/sigmax, StdErrRatioNormME.rat1, type = "l", xlim=c(0,0.3),                  
     xlab="", ylab="",main="",col=col.rat1, ylim =c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaNormME/sigmax, StdErrRatioNormME.rat2, type = "l", col = col.rat2, lty = lty.rat2)
lines(sddeltaNormME/sigmax, StdErrRatioNorm_n7_ME_term1only.rat3, type = "l",
      col = col.rat3, lty = lty.rat3)

text(pos.x.leg, pos.y.leg, expression(bold((b))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #





# Add an inner plot - Normal ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaNormME/sigmax, StdErrPercentNorm_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))
lines(sddeltaNormME/sigmax, StdErrPercentNorm_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
lines(sddeltaNormME/sigmax, StdErrPercentNorm_n7_ME.rat3, type = "l",lty =lty.rat3,
      xlab="",ylab="",col=col.rat3)
par(opar)
# end of 'add an inner plot - Normal ME'



# c. Beta measurement error ~~~~~~~~~~~
plot(sddeltaBetaME/sigmax, StdErrRatioBetaME.rat1, type ="l", xlim=c(0,0.3),                   
     xlab="", ylab="",  main="", col=col.rat1, ylim = c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaBetaME/sigmax, StdErrRatioBetaME.rat2, type = "l", col=col.rat2, lty = lty.rat2)
lines(sddeltaBetaME/sigmax, StdErrRatioBeta_n7_ME_term1only.rat3,
      type = "l", col=col.rat3, lty = lty.rat3)
text(pos.x.leg, pos.y.leg, expression(bold((c))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #

# Add an inner plot - Beta ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaBetaME/sigmax, StdErrPercentBeta_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)

lines(sddeltaBetaME/sigmax, StdErrPercentBeta_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
lines(sddeltaBetaME/sigmax, StdErrPercentBeta_n7_ME.rat3, type = "l",lty =lty.rat3,
      xlab="",ylab="",col=col.rat3)
par(opar)
# end of 'add an inner plot - Beta ME'
# end of 'Fig.8.2 - Std.Err' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# Fig.3.3 - RMSE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
y.lim.min=1;y.lim.max=2
y.lim.min.inn=3;y.lim.max.inn=7
loc.line.lab.y=2.2
pos.x.leg=0.29;pos.y.leg=y.lim.max*0.99;
pos.x.leg1=0.28;pos.y.leg1=y.lim.max*0.95;
lim.low.inn=8.7;lim.left.inn=mar.2;lim.up.inn=0.25;lim.rig.inn=6
per.pos.x.ylabel.inn=0.83;per.pos.y.ylabel.inn=0.98



# a. Uniform measurement error ~~~~~~~~~~~
plot(sddeltaUnifME/sigmax, RMSERatioUnif_n7_ME.rat1,col=col.rat1,
     type = "l", lty=1, xlim=c(0,0.3), xlab="",ylab="",main="",
     ylim = c(y.lim.min,y.lim.max)) 
lines(sddeltaUnifME/sigmax, RMSERatioUnif_n7_ME.rat2, type = "l", lty=lty.rat2,col=col.rat2)    
lines(sddeltaUnifME/sigmax, RMSERatioUnif_n7_ME.rat3, type = "l", lty=lty.rat3,col=col.rat3) 

title(ylab=list(expression("Raz\u00f3n de RECM"),col="brown"),
      line = loc.line.lab.y, font=2, cex.lab=1.0)
text(pos.x.leg, pos.y.leg, expression(bold((a))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #

# Add an inner plot - Uniform ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaUnifME/sigmax, RMSEpercentUnif_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1, ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))

lines(sddeltaUnifME/sigmax, RMSEpercentUnif_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
lines(sddeltaUnifME/sigmax, RMSEpercentUnif_n7_ME.rat3, type = "l",lty =lty.rat3,
      xlab="",ylab="",col=col.rat3)
par(opar)
# end of 'add an inner plot - Uniform ME'                             


# b. Normal measurement error ~~~~~~~~~~~
plot(sddeltaNormME/sigmax, RMSERatioNorm_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab="", ylab="",main="",col=col.rat1, ylim =c(y.lim.min, y.lim.max), lty=lty.rat1)
lines(sddeltaNormME/sigmax, RMSERatioNorm_n7_ME.rat2,type = "l", col=col.rat2, lty = lty.rat2)
lines(sddeltaNormME/sigmax,RMSERatioNorm_n7_ME.rat3, type = "l", col=col.rat3, lty = lty.rat3)

text(pos.x.leg, pos.y.leg, expression(bold((b))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #

# Add an inner plot - Normal ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaNormME/sigmax, RMSEpercentNorm_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))
lines(sddeltaNormME/sigmax, RMSEpercentNorm_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
lines(sddeltaNormME/sigmax, RMSEpercentNorm_n7_ME.rat3, type = "l",lty =lty.rat3,
      xlab="",ylab="",col=col.rat3)
par(opar)
# end of 'add an inner plot - Normal ME'



# c. Beta measurement error ~~~~~~~~~~~
plot(sddeltaBetaME/sigmax, RMSERatioBeta_n7_ME.rat1, type ="l", xlim=c(0,0.3),  
     xlab="", ylab="",  main="", col=col.rat1, ylim = c(y.lim.min, y.lim.max), lty=lty.rat1)
#abline(h= c(1,2), col="black", lty = 3)
lines(sddeltaBetaME/sigmax,RMSERatioBeta_n7_ME.rat2, type = "l", col=col.rat2, lty = lty.rat2)
lines(sddeltaBetaME/sigmax,RMSERatioBeta_n7_ME.rat3, type = "l", col=col.rat3, lty = lty.rat3)
text(pos.x.leg, pos.y.leg, expression(bold((c))))
text(pos.x.leg1*per.pos.x.ylabel.inn, pos.y.leg1*per.pos.y.ylabel.inn, "(%)") #

# Add an inner plot - Beta ME
opar <- par(new=T,mar=c(lim.low.inn,lim.left.inn,lim.up.inn,lim.rig.inn)) 
plot(sddeltaBetaME/sigmax, RMSEpercentBeta_n7_ME.rat1, type = "l", xlim=c(0,0.3),
     xlab=NA,ylab="", main="", col=col.rat1,
     ylim = c(y.lim.min.inn, y.lim.max.inn),
     lty=lty.rat1,cex.axis=1,axes =FALSE,frame.plot = TRUE)#xaxt="FALSE")
ys.used=seq(y.lim.min.inn,y.lim.max.inn,by=1)
axis(4, at = ys.used, cex.axis=0.8)#, labels = formatC(ys.used, format="fg"))
lines(sddeltaBetaME/sigmax, RMSEpercentBeta_n7_ME.rat2, type = "l",lty =lty.rat2,
      xlab="",ylab="",col=col.rat2)
lines(sddeltaBetaME/sigmax, RMSEpercentBeta_n7_ME.rat3, type = "l",lty =lty.rat3,
      xlab="",ylab="",col=col.rat3)
par(opar)
# end of 'add an inner plot - Beta ME'
# end of 'Fig.8.3 - RMSE' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# the outer botton label -----------
label.x.rand2=
  expression(paste("                    ",sep="   ", sigma [delta], " / ", sigma [x],
                   sep="                                         ",
                   sigma [delta], " / ", sigma [x],
                   sep="                                        ",
                   sigma [delta], " / ", sigma [x],
  ))

mtext(label.x.rand2,font=2,side=1,outer=TRUE,adj=0,line=0.9)

# end of the ouer botton label -------

par(oldpar)

Bibliografía

Gregoire, T. G., and C. Salas. 2009. “Ratio Estimation with Measurement Error in the Auxiliary Variate.” Biometrics 65 (2): 590–98.
Salas-Eljatib, C. 2021. Análisis de datos con el programa estadístico R: una introducción aplicada. Santiago, Chile: Ediciones Universidad Mayor.