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).

Última actualización de esta web: 13 noviembre, 2022

Sobre esta figura

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)

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)              
> #segregando los datos para facilitar los graficos luego
> 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")
> #definiendo las respectivas variables respuestas
> y   <- db90
> y70 <- db70
> y50 <- db50
> y30 <- db30
> y10 <- db10
> #definiendo los rangos para b0
> 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)
> #definiendo los rangos para b1
> 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)
> #definiendo los rangos para b2
> 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)
> #definiendo los rangos para b3
> 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)
> #definiendo los rangos para b4
> 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)
> #definiendo los rangos para b5
> 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)
> #reasignando y definiendo nuevos objetos
> 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. https://tienda.zigzag.cl/9789566086109-analisis-de-datos-con-el-programa-estadistico-r.html.
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. https://doi.org/10.1016/j.ecolind.2017.10.030.