Este código de R corresponde a la figura 6.16yy 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 30 histogramas con las siguientes características: (i) se marca en cada histograma una línea vertical; (ii) la etiqueta del eje-Y tiene dos renglones; (iii) la etiqueta del eje-X utiliza símbolos estadísticos con letras griegas.

Cargando los datos de ejemplo

La investigación de Salas-Eljatib et al. (2018) involucró estudiar la distribución de estimadores de máxima verosimilitud de modelos estadísticos de regresión logística. Esto se logro mediante simulaciones estocasticas. Los resultados se encuentran en el objeto de datos unbalanLogReg.RData que contienen valores estimados de todos los parámetros de un modelo estadístico de regresión logísticia.

library(datana)
load("unbalanLogReg.RData")
head(db)
  combo beta0.hat  beta1.hat  beta2.hat  beta3.hat beta4.hat  beta5.hat aic.out
1    90   -1.6644 -0.0041112 -0.0143511 -0.0026061   0.66893 0.00027465  617.53
2    90   -3.2371 -0.0023504  0.0875868 -0.0028754   0.85233 0.00014939  616.47
3    90   -1.8044 -0.0043650 -0.0079855 -0.0022202   0.84225 0.00026211  608.13
4    90   -1.9258 -0.0033947 -0.0804311 -0.0014277   0.80853 0.00030381  615.70
5    90   -3.6876 -0.0030000  0.0254684 -0.0033314   1.66340 0.00029861  584.88
6    90   -1.6358 -0.0054772  0.0487885 -0.0035931   1.00984 0.00026721  593.35
  deadBien deadMal NoDeadBien NoDeadMal Pred.correctas Pred.incorrectas
1        0     100        900         0             90               10
2        0     100        900         0             90               10
3        0     100        900         0             90               10
4        0     100        900         0             90               10
5        0     100        898         2             90               10
6        0     100        900         0             90               10

El gráfico

Ahora procedemos a realizar el gráfico.

#primero: un par de definiciones
col.bar <- "lightgray"
col.param <- "black"; lty.param <- 1 ;
lwd.param <- 2
col.exp.val <- "black"; lty.exp.val <- 2;
lwd.exp.val <- 2

num.param <- 6
ylab.freq <- "Frecuencia relativa"
num.breaks <- 20
ylim.h <- c(0, 0.2)


#segundo: la produccion de la figura
my.par <- par(mfrow = c(5,num.param),
 mar=c(0.0,3,1.5,0),oma=c(4,2,0,1),xaxs="i",yaxs="i",bty="l",
 cex.axis=0.6)              


db90 <- subset(db,combo=="90")
db70 <- subset(db,combo=="70")
db50 <- subset(db,combo=="50")
db30 <- subset(db,combo=="30")
db10 <- subset(db,combo=="10")

y   <- db90
y70 <- db70
y50 <- db50
y30 <- db30
y10 <- db10

low.b0 <- range(c(y$beta0.hat,y70$beta0.hat,y50$beta0.hat,y30$beta0.hat,y10$beta0.hat))[1]
hig.b0 <- range(c(y$beta0.hat,y70$beta0.hat,y50$beta0.hat,y30$beta0.hat,y10$beta0.hat))[2] 
breaks.b0 <- seq(low.b0, hig.b0, length=num.breaks)


low.b1 <- range(y$beta1.hat)[1]
hig.b1 <- range(y$beta1.hat)[2]
low.b1 <- range(c(y$beta1.hat,y70$beta1.hat,y50$beta1.hat,y30$beta1.hat,y10$beta1.hat))[1]
hig.b1 <- range(c(y$beta1.hat,y70$beta1.hat,y50$beta1.hat,y30$beta1.hat,y10$beta1.hat))[2]
breaks.b1 <- seq(low.b1, hig.b1, length=num.breaks)

low.b2 <-  range(c(y$beta2.hat,y70$beta2.hat,y50$beta2.hat,y30$beta2.hat,y10$beta2.hat))[1]
hig.b2 <-  range(c(y$beta2.hat,y70$beta2.hat,y50$beta2.hat,y30$beta2.hat,y10$beta2.hat))[2]    
breaks.b2 <- seq(low.b2, hig.b2, length=num.breaks)

low.b3 <- range(c(y$beta3.hat,y70$beta3.hat,y50$beta3.hat,y30$beta3.hat,y10$beta3.hat))[1]
hig.b3 <- range(c(y$beta3.hat,y70$beta3.hat,y50$beta3.hat,y30$beta3.hat,y10$beta3.hat))[2]
breaks.b3 <- seq(low.b3, hig.b3, length=num.breaks)

low.b4 <- range(c(y$beta4.hat,y70$beta4.hat,y50$beta4.hat,y30$beta4.hat,y10$beta4.hat))[1]
hig.b4 <- range(c(y$beta4.hat,y70$beta4.hat,y50$beta4.hat,y30$beta4.hat,y10$beta4.hat))[2]
breaks.b4 <- seq(low.b4, hig.b4, length=num.breaks)

low.b5 <-range(c(y$beta5.hat,y70$beta5.hat,y50$beta5.hat,y30$beta5.hat,y10$beta5.hat))[1]
hig.b5 <- range(c(y$beta5.hat,y70$beta5.hat,y50$beta5.hat,y30$beta5.hat,y10$beta5.hat))[2]
breaks.b5 <- seq(low.b5, hig.b5, length=num.breaks)

valor.minx.b0 <-  low.b0; valor.maxx.b0 <-  hig.b0
valor.minx.b1 <-  low.b1; valor.maxx.b1 <-  hig.b1
valor.minx.b2 <-  low.b2; valor.maxx.b2 <-  hig.b2
valor.minx.b3 <-  low.b3; valor.maxx.b3 <-  hig.b3
valor.minx.b4 <-  low.b4; valor.maxx.b4 <-  hig.b4
valor.minx.b5 <-  low.b5; valor.maxx.b5 <-  hig.b5

#Sexta linea tmo 90%
main.h='90% de ceros'
y   <- db90

#**************
#B0-90
#**************
h <- hist(y$beta0.hat, plot=F, breaks=breaks.b0)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b0, valor.maxx.b0), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",xaxt='n', las=1)
abline(h=0)

abline(v=b0.real, col=col.param,lty=lty.param, lwd=lwd.param )

b.param.h <- b0.real
y.here <- y$beta0.hat
exp.h <- mean(y.here)

exp.b0.90 <- exp.h

abline(v=exp.h, col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)
mtext(main.h,side=2,line=3,font=1,cex=0.8)
mtext(ylab.freq,side=2,line=2.2,font=1,cex=0.5)

#esto es para la segunda columna
par(mar=c(0.0,0.9,1.5,0))

#**************
#B1-90
#**************
h = hist(y$beta1.hat, plot=F, breaks=breaks.b1)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b1, valor.maxx.b1), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b1.real, col=col.param,lty=lty.param, lwd=lwd.param )


b.param.h <- b1.real
y.here <- y$beta1.hat
exp.h <- mean(y.here)

exp.b1.90 <- exp.h

abline(v=exp.h, col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

#**************
##B2-90
#**************
h = hist(y$beta2.hat, plot=F, breaks=breaks.b2)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b2, valor.maxx.b2), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

b.param.h <- b2.real
y.here <- y$beta2.hat
exp.h <- mean(y.here)

exp.b2.90 <- exp.h

abline(v=b2.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta2.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

#**************
#B3-90
#**************
h = hist(y$beta3.hat, plot=F, breaks=breaks.b3)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b3, valor.maxx.b3), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b3.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta3.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b3.real
y.here <- y$beta3.hat
exp.h <- mean(y.here)

exp.b3.90 <- exp.h

#**************
#B4-90
#**************
h = hist(y$beta4.hat, plot=F, breaks=breaks.b4)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="",
     xlim=c(valor.minx.b4, valor.maxx.b4),
     ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b4.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta4.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b4.real
y.here <- y$beta4.hat
exp.h <- mean(y.here)

exp.b4.90 <- exp.h

#**************
#B5-90
#**************
h = hist(y$beta5.hat, plot=F, breaks=breaks.b5)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b5, valor.maxx.b5), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b5.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta5.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b5.real
y.here <- y$beta5.hat
exp.h <- mean(y.here)

exp.b5.90 <- exp.h

par(mar=c(0.0,3,1.5,0))

#=======================
#Segunda linea de tmo 70%
message("---------------------------")
---------------------------
message("Plotting for 70% of zeros")
Plotting for 70% of zeros
y <- db70
main.h='70% de ceros'

#**************
#B0-70
#**************
h = hist(y$beta0.hat, plot=F, breaks=breaks.b0)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b0, valor.maxx.b0), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",xaxt='n',las=1)
abline(h=0)

abline(v=b0.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta0.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b0.real
y.here <- y$beta0.hat
exp.h <- mean(y.here)

exp.b0.70 <- exp.h

mtext(main.h,side=2,line=3,font=1,cex=0.8)
mtext(ylab.freq,side=2,line=2.2,font=1,cex=0.5)

par(mar=c(0.0,0.9,1.5,0))
#@@@@@@@@
#B1-70
#@@@@@@@@
h = hist(y$beta1.hat, plot=F, breaks=breaks.b1)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b1, valor.maxx.b1), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

b.param.h <- b1.real
y.here <- y$beta1.hat
exp.h <- mean(y.here)

exp.b1.70 <- exp.h

abline(v=b1.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta1.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

##B2-70
h = hist(y$beta2.hat, plot=F, breaks=breaks.b2)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b2, valor.maxx.b2), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

b.param.h <- b2.real
y.here <- y$beta2.hat
exp.h <- mean(y.here)

exp.b2.70 <- exp.h

abline(v=b2.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta2.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

#B3-70
h = hist(y$beta3.hat, plot=F, breaks=breaks.b3)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b3, valor.maxx.b3), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b3.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta3.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b3.real
y.here <- y$beta3.hat
exp.h <- mean(y.here)

exp.b3.70 <- exp.h

#B4-70
h = hist(y$beta4.hat, plot=F, breaks=breaks.b4)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b4, valor.maxx.b4), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b4.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta4.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b4.real
y.here <- y$beta4.hat
exp.h <- mean(y.here)

exp.b4.70 <- exp.h

#**************
#B5-70
#**************
h = hist(y$beta5.hat, plot=F, breaks=breaks.b5)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b5, valor.maxx.b5), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b5.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta5.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b5.real
y.here <- y$beta5.hat
exp.h <- mean(y.here)

exp.b5.70 <- exp.h

#esto es para la segunda linea, vuelvo a los mar() originales
par(mar=c(0.0,3,1.5,0))

#======================
#Tercera linea tmo 50%
message("---------------------------")
---------------------------
message("Plotting for 50% of zeros")
Plotting for 50% of zeros
y <- db50

main.h='50% de ceros'

#B0-50
h = hist(y$beta0.hat, plot=F, breaks=breaks.b0)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b0, valor.maxx.b0), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",xaxt='n',las=1)
abline(h=0)

abline(v=b0.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta0.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b0.real
y.here <- y$beta0.hat
exp.h <- mean(y.here)

exp.b0.50 <- exp.h

mtext(main.h,side=2,line=3,font=1,cex=0.8)
mtext(ylab.freq,side=2,line=2.2,font=1,cex=0.5)


par(mar=c(0.0,0.9,1.5,0))
#B1-50
h = hist(y$beta1.hat, plot=F, breaks=breaks.b1)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b1, valor.maxx.b1), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",xaxt='n',yaxt='n')
abline(h=0)

abline(v=b1.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta1.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b1.real
y.here <- y$beta1.hat
exp.h <- mean(y.here)

exp.b1.50 <- exp.h

##B2-50
h = hist(y$beta2.hat, plot=F, breaks=breaks.b2)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b2, valor.maxx.b2), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b2.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta2.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b2.real
y.here <- y$beta2.hat
exp.h <- mean(y.here)

exp.b2.50 <- exp.h

#B3-50
h = hist(y$beta3.hat, plot=F, breaks=breaks.b3)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b3, valor.maxx.b3), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)


b.param.h <- b3.real
y.here <- y$beta3.hat
exp.h <- mean(y.here)

exp.b3.50 <- exp.h

abline(v=b3.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta3.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

#B4-50
h = hist(y$beta4.hat, plot=F, breaks=breaks.b4)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b4, valor.maxx.b4),
     ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b4.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta4.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b4.real
y.here <- y$beta4.hat
exp.h <- mean(y.here)

exp.b4.50 <- exp.h

#**************
#B5-50
#**************
h = hist(y$beta5.hat, plot=F, breaks=breaks.b5)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b5, valor.maxx.b5), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b5.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta5.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b5.real
y.here <- y$beta5.hat
exp.h <- mean(y.here)

exp.b5.50 <- exp.h

#esto es para la segunda linea, vuelvo a los mar() originales
par(mar=c(0.0,3,1.5,0))
#==========================
# Pen-Ultima linea de tmo 30%
#Segunda linea tmo 30%
y <- db30

message("---------------------------")
---------------------------
message("Plotting for 30% of zeros")
Plotting for 30% of zeros
main.h='30% de ceros'

#B0-30
h = hist(y$beta0.hat, plot=F, breaks=breaks.b0)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b0, valor.maxx.b0), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",xaxt='n',las=1)
abline(h=0)
abline(v=b0.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta0.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b0.real
y.here <- y$beta0.hat
exp.h <- mean(y.here)

exp.b0.30 <- exp.h

mtext(main.h,side=2,line=3,font=1,cex=0.8)
mtext(ylab.freq,side=2,line=2.2,font=1,cex=0.5)

#esto es para la segunda columna
par(mar=c(0.0,0.9,1.5,0))
#B1-30
h = hist(y$beta1.hat, plot=F, breaks=breaks.b1)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b1, valor.maxx.b1), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b1.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta1.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b1.real
y.here <- y$beta1.hat
exp.h <- mean(y.here)
exp.b1.30 <- exp.h

##B2-30
h = hist(y$beta2.hat, plot=F, breaks=breaks.b2)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b2, valor.maxx.b2), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b2.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta2.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b2.real
y.here <- y$beta2.hat
exp.h <- mean(y.here)

exp.b2.30 <- exp.h

#B3-30
h = hist(y$beta3.hat, plot=F, breaks=breaks.b3)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b3, valor.maxx.b3), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b3.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta3.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b3.real
y.here <- y$beta3.hat
exp.h <- mean(y.here)

exp.b3.30 <- exp.h

#B4-30
h = hist(y$beta4.hat, plot=F, breaks=breaks.b4)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b4, valor.maxx.b4), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)
abline(v=b4.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta4.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b4.real
y.here <- y$beta4.hat
exp.h <- mean(y.here)

exp.b4.30 <- exp.h


#**************
#B5-30
#**************
h = hist(y$beta5.hat, plot=F, breaks=breaks.b5)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b5, valor.maxx.b5), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n',xaxt='n')
abline(h=0)

abline(v=b5.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta5.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b5.real
y.here <- y$beta5.hat
exp.h <- mean(y.here)

exp.b5.30 <- exp.h

#esto es para la segunda linea, vuelvo a los mar() originales
par(mar=c(0.0,3,1.5,0))
#==========================
# Ultima linea de tmo 10%
y <- db10

message("---------------------------")
---------------------------
message("Plotting for 10% of zeros")
Plotting for 10% of zeros
main.h='10% de ceros'

#B0-10
h = hist(y$beta0.hat, plot=F, breaks=breaks.b0)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b0, valor.maxx.b0), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",las=1)
abline(h=0)
abline(v=b0.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta0.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b0.real
y.here <- y$beta0.hat
exp.h <- mean(y.here)

exp.b0.10 <- exp.h

mtext(main.h,side=2,line=3,font=1,cex=0.8)
mtext(ylab.freq,side=2,line=2.2,font=1,cex=0.5)
mtext(expression(widehat(beta)[0]),side=1,line=3,font=1, cex=1)

#esto es para la segunda columna
par(mar=c(0.0,0.9,1.5,0))
#B1-10
h = hist(y$beta1.hat, plot=F, breaks=breaks.b1)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b1, valor.maxx.b1), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n')
abline(h=0)
abline(v=b1.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta1.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b1.real
y.here <- y$beta1.hat
exp.h <- mean(y.here)

exp.b1.10 <- exp.h

mtext(expression(widehat(beta)[1]),side=1,line=3,font=1, cex=1)

##B2-10
h = hist(y$beta2.hat, plot=F, breaks=breaks.b2)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b2, valor.maxx.b2), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n')
abline(h=0)
abline(v=b2.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta2.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b2.real
y.here <- y$beta2.hat
exp.h <- mean(y.here)

exp.b2.10 <- exp.h

mtext(expression(widehat(beta)[2]),side=1,line=3,font=1, cex=1)


#B3-10
h = hist(y$beta3.hat, plot=F, breaks=breaks.b3)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="", xlim=c(valor.minx.b3, valor.maxx.b3), ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n')
abline(h=0)
abline(v=b3.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta3.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b3.real
y.here <- y$beta3.hat
exp.h <- mean(y.here)

exp.b3.10 <- exp.h

mtext(expression(widehat(beta)[3]),side=1,line=3,font=1, cex=1) 

#B4-10
h = hist(y$beta4.hat, plot=F, breaks=breaks.b4)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="",
     xlim=c(valor.minx.b4, valor.maxx.b4),
     ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n')
abline(h=0)
abline(v=b4.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta4.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b4.real
y.here <- y$beta4.hat
exp.h <- mean(y.here)

exp.b4.10 <- exp.h

mtext(expression(widehat(beta)[4]),side=1,line=3,font=1, cex=1)


#B5-10
h = hist(y$beta5.hat, plot=F, breaks=breaks.b5)
h$density <- h$density/sum(h$density) 
plot(h, freq=F, xlab="", ylab="",
     xlim=c(valor.minx.b5, valor.maxx.b5),
     ylim=ylim.h, border=FALSE, col=col.bar,
     main="",yaxt='n')
abline(h=0)
abline(v=b5.real, col=col.param,lty=lty.param, lwd=lwd.param )
abline(v=mean(y$beta5.hat), col=col.exp.val,lty=lty.exp.val,
       lwd=lwd.exp.val)

b.param.h <- b5.real
y.here <- y$beta5.hat
exp.h <- mean(y.here)

exp.b5.10 <- exp.h

mtext(expression(widehat(beta)[5]),side=1,line=3,font=1, cex=1)

Bibliografía

Salas-Eljatib, C. 2021. Análisis de datos con el programa estadístico R: una introducción aplicada. Santiago, Chile: Ediciones Universidad Mayor.
Salas-Eljatib, C., A. Fuentes-Ramírez, T. G. Gregoire, A. Altamirano, and V. Yaitul. 2018. “A Study on the Effects of Unbalanced Data When Fitting Logistic Regression Models in Ecology.” Ecological Indicators 85: 502–8.