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

Este tipo de gráfico mezcla no tan sólo observaciones sino que también la predicción de modelos e intervalos de confianza alrededor de dichas predicciones. Los datos y modelos provienen del estudio de Soto et al. (2015).

Cargando los datos de ejemplo

Se ocupará el objeto de datos modseedlingrauli.RData que es una lista conteniendo datos y vectores necesarios para realizar el gráfico. Este objeto se encuentra disponible en el paquete datana. Los datos y modelos provienen del estudio de Soto et al. (2015).

library(datana)
load("modseedlingrauli.RData")
head(db0)#revisemos una de las dataframes importantes
    LTR     GTR   LIR     GIR   LBR     GBR   LAR     GAR   LTC     GTC   LIC
1 18.35 0.55962 41.16 0.36772 18.35 0.55962 14.21 0.10536 44.20 0.61904 51.02
2 33.75 0.59784 32.60 0.55962 33.75 0.59784 21.53 0.33647 27.24 0.40547 54.33
3 43.19 0.89382 36.53 0.62861 43.19 0.89382 67.17 0.53063 36.47 0.53900 25.84
4 44.15 0.62861 27.80 0.57536 44.15 0.62861 70.91 0.45861 40.58 0.77319 26.21
5 40.95 0.89382 36.31 0.84730 40.95 0.89382 19.54 0.26236 24.45 0.40547 29.81
6 18.59 0.40547 41.56 0.84730 18.59 0.40547 45.68 0.40547 31.05 0.69315 41.49
      GIC   LBC     GBC   LAC     GAC    y.co    y.ra  x.co  x.ra
1 0.53900 44.20 0.61904 40.65 0.60614 0.61904 0.55962 44.20 18.35
2 0.69315 27.24 0.40547 33.29 0.74194 0.40547 0.59784 27.24 33.75
3 0.91629 36.47 0.53900 40.13 0.76214 0.53900 0.89382 36.47 43.19
4 0.61904 40.58 0.77319 47.25 0.40547 0.77319 0.62861 40.58 44.15
5 0.47000 24.45 0.40547 67.28 0.57536 0.40547 0.89382 24.45 40.95
6 0.55962 31.05 0.69315 40.39 0.28768 0.69315 0.40547 31.05 18.59

El gráfico

Ahora se procede a realizar el gráfico, pero se explica por etapas para así comprenderlo de mejor forma.

Primero se define los nombres, colores, y tipos de líneas, a ser usados en la leyenda del gráfico

pch.l<-3;pch.i<-4;pch.h<-10
name.l<-"Baja";name.i<-"Media";name.h<-"Alta"
name.lm<-"Mod. comp. baja";name.im<-"Mod. comp. media";
name.hm<-"Mod. comp. alta"
col.l<-"blue";col.i<-"black";col.h<-"red"
col.envl<-"lightblue";col.envi<-"lightgray";col.envh<-"pink"
lty.l<-2;lty.i<-1;lty.h<-4
lwd.l=1;lwd.i=1;lwd.h=2;

legend.names <- c(name.l,name.i,name.h,name.lm,name.im,name.hm)#,name.mod3)
my.col.names <- c(col.l,col.i,col.h,col.l,col.i,col.h)

Un par de objetos de importancia también son

x.list <- seq(from=0, to=80, by=.1)

head(pred)
        x
1 0.00000
2 0.80808
3 1.61616
4 2.42424
5 3.23232
6 4.04040
head(pred2.ra.h)
    y.pred     low.y    upp.y
1 0.000000 0.0000000 0.000000
2 0.014850 0.0058111 0.023888
3 0.029155 0.0120056 0.046304
4 0.042946 0.0185299 0.067362
5 0.056249 0.0253365 0.087162
6 0.069090 0.0323832 0.105797
head(data.frame(d.ra.l,d.ra.i,d.ra.h))
     d.ra.l    d.ra.i    d.ra.h
1 0.0000000 0.0000000 0.0000000
2 0.0027552 0.0027691 0.0018682
3 0.0055023 0.0055283 0.0037276
4 0.0082412 0.0082778 0.0055784
5 0.0109721 0.0110176 0.0074205
6 0.0136949 0.0137476 0.0092540
par(xaxs="i", yaxs="i", bty='l')#,col.axis='white')
plot(GBR~LBR, data=db0, las=1, col=col.l, 
     xlab="",  
     ylab="",pch=0,cex=0,
     xaxt="n", yaxt="n", 
     ylim=c(0,1.4),xlim=c(0,80))

#CI envelope for Ra.high
polygon(c(pred$x,rev(pred$x)),c(pred2.ra.h[,2],rev(pred2.ra.h[,3])),
        col=col.envh,border=NA) #

#CI envelope for Ra.low
polygon(c(pred$x,rev(pred$x)),c(pred2.ra.l[,2],rev(pred2.ra.l[,3])),
        col=col.envl,border=NA) #

#CI envelope for Ra.int
polygon(c(pred$x,rev(pred$x)),c(pred2.ra.i[,2],rev(pred2.ra.i[,3])),
        col=col.envi,border=NA) #

par(new=T)
plot(GBR~LBR, data=db0, las=1, col=col.l, 
     xlab="Radiaci\u00f3n transmitida (%)",  
     ylab="Incremento relativo en di\u00E1metro",pch=pch.l, cex=0.6,
     ylim=c(0,1.4),xlim=c(0,80))
points(db0$LIR, db0$GIR,  cex=0.6, col=col.i, pch=pch.i)
points(db0$LAR, db0$GAR,  cex=0.6, col=col.h, pch=pch.h)

#la prediccion, o valor esperado, de cada modelo
lines(d.ra.l~x.list,col=col.l,lty=lty.l,lwd=lwd.l)
lines(d.ra.i~x.list,col=col.i,lty=lty.i,lwd=lwd.i)
lines(d.ra.h~x.list,col=col.h,lty=lty.h,lwd=lwd.h)

legend("topleft",legend.names,bty="n",
       lty=c(NA,NA,NA,lty.l,lty.i,lty.h),
       lwd=c(NA,NA,NA,lwd.l,lwd.i,lwd.h),
       pch=c(pch.l,pch.i,pch.h,NA,NA,NA),
       col=my.col.names)