Mapear archivo KML

Con datos de la Secretaría de Energía (SENER) se mapeará una capa de información en formato KML (del acrónimo en inglés Keyhole Markup Language). Los datos corresponden a las zonas en proceso de electrificación.

El código para extraer las coordenadas geográficas, la altitud así como el identificador (ID), es un aporte de Fabio Veronesi (https://www.r-bloggers.com/extract-coordinates-and-other-data-from-kml-in-r/), ligeras modificaciones fueron efectuadas. No obstante como veremos, los datos de la SENER sólo contienen latitud y longitud, la forma de hacerlo es la siguiente.

kml.text <- readLines('SENER.kml')
re <- "<coordinates> *([^<]+?) *<\\/coordinates>"  
coords <- grep(re,kml.text)  
   
re2 <- "src_id:"  
SCR.ID <- grep(re2,kml.text)  
   
re3 <- "<tr><td><b>Name:</b><td>"  
Name <- grep(re3,kml.text)  
   
kml.coordinates <- matrix(0,length(coords),4,dimnames=list(c(),c("ID","LON","LAT","ELEV")))  
kml.names <- matrix(0,length(coords),1)  
   
for(i in 1:length(coords)){  
    sub.coords <- coords[i]  
    temp1 <- gsub("<coordinates>"," ",kml.text[sub.coords])  
    temp2 <- gsub("</coordinates>"," ",temp1)  
    coordinates <- as.numeric(unlist(strsplit(temp2,",")))  
   
    sub.ID <- SCR.ID[i]  
    ID <- as.numeric(gsub("<tr><td><b>src_id:</b><td>"," ",kml.text[sub.ID]))  
   
    sub.Name <- Name[i]  
    NAME <- gsub(paste("<tr><td><b>Name:</b><td>"),"",kml.text[sub.Name])  
   
    kml.coordinates[i,] <- matrix(c(ID,coordinates),ncol=4)  
    kml.names[i,] <- matrix(c(NAME),ncol=1)  
}

Con esta parte del código extraemos la información del archivo KML.

> head(kml.coordinates)
     ID       LON      LAT ELEV
[1,] NA -92.13972 17.18167    0
[2,] NA -91.99667 17.18222    0
[3,] NA -91.20667 16.37694    0
[4,] NA -91.74639 17.03611    0
[5,] NA -92.11972 17.09611    0
[6,] NA -92.80972 17.40472    0

Como vemos, la información proporcionada por la dependencia se limita a las coordenadas geográficas. Por lo que, en este ejemplo, dejaremos la matriz sólo con las columnas que contienen información.

kml.coordinates <- kml.coordinates[,2:3]

El siguiente paso será plotear estas coordenadas, el código empleado es el siguiente:

library(RgoogleMaps)
z0 <- 5                     # zoom
y0 <- 21.525818             # Latitud zero
x0 <- -100.522399           # longitud zero


N <- length(kml.coordinates)
color <- c()                # Inicializo vector
color[1:N] <- 'gray20'

PuebMap  <- GetMap(center= c(y0,x0), zoom=z0, 
             GRAYSCALE = F, size=c(640,640),
             destfile = file.path(tempdir(),"sener.png"),
             maptype = 'terrain', SCALE = 2,)

PlotOnStaticMap(PuebMap)

{PlotOnStaticMap(PuebMap, lat = kml.coordinates[,2], 
                lon = kml.coordinates[,1], cex=1, pch=19,
                destfile = "sener.png", 
                col= color ,
                add=TRUE)};

Esto produce el siguiente mapa:

Esta representación nos permite observar que los procesos de electrificación se están llevando a cabo principalmente en las áreas montañosas.

Análisis Pokemon

La base de datos fue reunida de diversas fuentes por Alberto Barradas (https://www.kaggle.com/abcsds), la cual incluye 721 pokemon con algunas de sus estadísticas básicas:

> pokemon <- as.data.frame(read.csv('Pokemon.csv', header=T))
> str(pokemon)
'data.frame':   800 obs. of  13 variables:
 $ X.        : int  1 2 3 3 4 5 6 6 6 7 ...
 $ Name      : Factor w/ 800 levels "Abomasnow","AbomasnowMega Abomasnow",..: 81 330 746 747 103 104 100 101 102 666 ...
 $ Type.1    : Factor w/ 18 levels "Bug","Dark","Dragon",..: 10 10 10 10 7 7 7 7 7 18 ...
 $ Type.2    : Factor w/ 19 levels "","Bug","Dark",..: 15 15 15 15 1 1 9 4 9 1 ...
 $ Total     : int  318 405 525 625 309 405 534 634 634 314 ...
 $ HP        : int  45 60 80 80 39 58 78 78 78 44 ...
 $ Attack    : int  49 62 82 100 52 64 84 130 104 48 ...
 $ Defense   : int  49 63 83 123 43 58 78 111 78 65 ...
 $ Sp..Atk   : int  65 80 100 122 60 80 109 130 159 50 ...
 $ Sp..Def   : int  65 80 100 120 50 65 85 85 115 64 ...
 $ Speed     : int  45 60 80 80 65 80 100 100 100 43 ...
 $ Generation: Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Legendary : Factor w/ 2 levels "False","True": 1 1 1 1 1 1 1 1 1 1 ...

• X.: ID para cada pokemon
• Name: Nombre de cada pokemon
• Type.1: Cada pokemon tiene un tipo, el cual determina vulnerabilidades o fortalezas a los ataques.
• Type.2; Algunos pokemones son tipo dual y tienen un segundo "tipo".
• Total: Suma de las estadísticas que se listan a continuación, da una idea general sobre que tan fuerte es el pokemon.
• HP: Hit points, o salud, define que tanto daño puede soportar un pokemon después de pelear
• Attack: Modificador de base para ataques normales.
• Defense: Resistencia al daño base contra ataques normales.
• SP Atk: Ataques especiales, el modificador de base para ataques especiales como rafaga de fuego, etc.
• SP Def: Resistencia al daño base contra ataques especiales.
• Speed: Determina que pokemon ataca primero cada ronda.
• Generation: La generación a la que pertenece el pokemon.
• Legendary: Variable booleana que indica TRUE cuando el pokemon es catalogado como legendario.

Promedio "Total" por tipo (1) de pokemon y generación

pkmn <- melt(with(pokemon, tapply(Total, list(Generation,Type.1), mean)))

pkmn <- pkmn[,c(2,1,3)]

ina  <- which(is.na(pkmn$value))

pkmn <- pkmn[-ina,]

colnames(pkmn) <- c('Type','Generation','Total')

freq <- count(pokemon, c('Type.1','Generation'))

pkmn$Freq <- freq[,3]

rownames(pkmn) <- seq(1,nrow(pkmn))

pkmn

p <- ggplot(pkmn) +

    geom_point(aes(x      = as.character(Type),

                   y      = as.character(Generation),

                   size   = Freq, 

                   color  = Total)) +

    scale_x_discrete(name = 'Type') +

    scale_y_discrete(name = 'Generation') +

    scale_colour_gradient(name='Avg Total', low='yellow', high='red3') +

    theme_bw() + 

    theme(axis.text.x=element_text(angle = 90, hjust = 0, vjust = 0.5)) + 

    ggtitle('Pokemons by Generation, type \n and Total Score')

p

Los pokemon tipo "Dragón" tienen en promedio los mejores puntajes, por otro lado los pokemones tipo "Agua" son los más abundantes.

Tipos de pokemon y puntaje "Total"

pkmn.1 <- melt(with(pokemon, tapply(Total, Type.1, mean )))

colnames(pkmn.1) <- c('Type.1','Total')

pkmn.1$freq <- melt(table(pokemon$Type.1))[,2]

mypalette <- rev(rainbow(20))

treemap(pkmn.1, index=c('Type.1'), vSize='freq', vColor='Total', 

        palette= mypalette,type='value', title.legend='Avg Total', 

        title='Pokemons by type and Total')

Clasificando pokemones legendarios

Convertimos la variable "Generation" a escala ordinal

pokemon$Generation <- as.ordered(pokemon$Generation)

Dividimos la base de datos original, en una muestra de entrenamiento con el 60% de los datos y el resto será la muestra de evaluación del modelo:

set.seed(1111)

id.train <- sample(seq(1,nrow(pokemon)), round(nrow(pokemon)*0.6))

train <- pokemon[id.train,]        # Muestra entrenamiento

test  <- pokemon[-id.train,]       # Muestra evaluación

El modelo sólo considerará las variables numéricas y la variable ordinal que definimos previamente.

> form <- as.formula(paste('Legendary', paste(names(train[c(5:12)]), 

+                    collapse=' + '), sep=' ~ '))

> form

Legendary ~ Total + HP + Attack + Defense + Sp..Atk + Sp..Def + 

    Speed + Generation

Regresión logística

> log <- glm(form, data=train, family = binomial())

> summary(log)

Call:

glm(formula = form, family = binomial(), data = train)

Deviance Residuals: 

     Min        1Q    Median        3Q       Max  

-2.41554  -0.14484  -0.03142  -0.00470   2.42037  

Coefficients: (1 not defined because of singularities)

               Estimate Std. Error z value Pr(>|z|)    

(Intercept)  -22.619056   3.338589  -6.775 1.24e-11 ***

Total          0.047109   0.011872   3.968 7.25e-05 ***

HP            -0.021168   0.015201  -1.393   0.1638    

Attack        -0.026050   0.015030  -1.733   0.0831 .  

Defense       -0.008286   0.011818  -0.701   0.4832    

Sp..Atk       -0.005889   0.014302  -0.412   0.6805    

Sp..Def       -0.002507   0.013974  -0.179   0.8576    

Speed                NA         NA      NA       NA    

Generation.L   1.660818   0.768322   2.162   0.0306 *  

Generation.Q  -0.438145   0.718330  -0.610   0.5419    

Generation.C   0.113477   0.711620   0.159   0.8733    

Generation^4  -0.715514   0.640430  -1.117   0.2639    

Generation^5  -0.565000   0.570013  -0.991   0.3216    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 284.85  on 479  degrees of freedom

Residual deviance: 112.57  on 468  degrees of freedom

AIC: 136.57

Number of Fisher Scoring iterations: 8

La fortaleza general del pokemon es el mejor predictor de los pokemones legendarios, seguido del ataque, y en menor medida la resistencia al daño. De igual forma, la primera generación de pokemones, se distingue significativamente de la 6a generación en cuanto a la cantidad de legendarios.

> cm.l <- table(train$Legendary, ifelse(fitted(log)>0.5,1,0));cm.l

       

          0   1

  False 430   8

  True   17  25

> cr.l <- sum(diag(cm.l)/sum(cm.l));cr.l             # Razón de exactitud

[1] 0.9479167

> rsq<- (1 - (log$deviance/log$null.deviance));rsq   #Pseudo R cuadrada

[1] 0.6047947

La matriz de confusión sobre la muestra de entrenamiento, y la razón de "exactitud" indican que es un buen modelo, no así el valor de la pseudo-R cuadrada. En cuanto al pronóstico sobre la muestra de evaluación:

> cm.l <- table(test$Legendary, ifelse(test$score.l>0.5,1,0));cm.l

       

          0   1

  False 292   5

  True   10  13

> cr.l <- sum(diag(cm.l)/sum(cm.l));cr.l

[1] 0.953125

>    pred.l      <- prediction(test$score.l,test$Legendary)

>    perf.l      <- performance(pred.l,"tpr","fpr")

>    max(attr(perf.l,'y.values')[[1]]-attr(perf.l,'x.values')[[1]])

[1] 0.9360269

> auc.l <- paste('Logistic=',round(performance(pred.l,'auc')@y.values[[1]],3))

> auc.l

[1] "Logistic= 0.984"

De acuerdo con la matriz de confusión, la razón de "exactitud", el estadístico de Kolmogórov-Smirnov y el Área Bajo la Curva, se puede apreciar que es un buen modelo.

A pesar de lo bueno que resulta ser el modelo de clasificación de máxima entropia, o de regresión logística, en el siguiente gráfico se muestran los resultados del desempeño predictivo de otros dos modelos de clasificación: Splines de Regresión Adaptativa Multivariante (MARS por sus siglas en inglés) y AdaBoost con el algoritmo SAMME.

clr <- palette(rev(rainbow(3)));clr <- palette(rev(rainbow(3)))

plot(perf.l, col=clr[1], lty=6, lwd=2,        # LOG

     xlab='Tasa Falso Positivo',ylab='Tasa Verdadero Positivo')

plot(perf.m, add=T, col=clr[2],lty=5, lwd=2)  # MARS

plot(perf.a, add=T, col=clr[3],lty=4, lwd=2)  # ADABOOST

abline(a=0, b= 1)

legend(0.20,0.80,                                 

      c(auc.l,auc.m,auc.a),                       

      col=clr, lwd=3, lty= c(6,5,4))              

grid(col='gray85')

GoogleMaps + Bubble plot en R

En esta entrada muestro cómo graficar referencias geográficas empleando GoogleMaps y R. Para empezar necesitamos referencias geográficas que deseemos plotear, en mi caso graficaré unidades económicas de la Comisión Federal de Electricidad de los estados de Tlaxcala y Puebla, los datos los obtuve del Directorio Estadístico Nacional de Unidades Económicas (DENUE) del INEGI.

Instalamos y cargamos la librería

install.packages('RgoogleMaps')
library(RgoogleMaps)

Posterior a ello, definimos algunos parámetros como el zoom de nuestro mapa, la latitud y la longitud centro, y adicional a ello una distancia (Euclidiana) que será el radio de círculos concéntricos a partir de nuestro centro.

z0 <- 8                     # zoom
y0 <- 19.040442             # Latitud zero
x0 <- -98.250160            # longitud zero
k  <- 0.4                   # Define el salto del circulo
                            # Saltos de 11 Km por cada décima

Colocaremos en el centro del mapa (o cualquier coordenada), un marcador que nos indique alguna referencia adicional.

markers = paste0("&markers=color:blue|label:I|19.040442,-98.250160")

Acto seguido, obtener nuestro mapa base:

PuebMap  <- GetMap(center= c(y0,x0), zoom=z0, GRAYSCALE = TRUE,
             size=c(640,640),markers=markers,
             destfile = file.path(tempdir(),"interm.png"),
             maptype = 'mapmaker-hybrid')

Enseguida, importamos nuestra base de datos de las unidades económicas, la cual deberá tener una estructura parecida a esto:

> head(DF)
       id      lat       lon        x1 x2
1 3535207 18.87717 -97.72655   87.5724  3
2 3398690 18.49201 -97.40435 1167.6320 40
3 3263793 19.04774 -98.04557   87.5724  3
4 3467483 19.12102 -98.24399 1167.6320 40
5 3316191 18.95277 -97.91986   87.5724  3
6 3457112 19.06409 -98.23185   87.5724  3

x1 y x2 representan variables cuantitativas a representar mediante un gráfico de burbujas. Lo siguiente es elaborar el gráfico sobre nuestro mapa base.

tpm <- bubbleMap(DF,coords=c("lon","lat"),

                 crs=sp::CRS("+proj=longlat +datum=WGS84"),

                 max.radius=500,legendLoc="bottomright",

                 filename = file.path(tempdir(),"interm.png"),

                 map=PuebMap,zcol='zinc',

                 key.entries = c(100,250,500,750,2500,6000));

Lucirá de la siguiente manera:

Lo siguiente que haremos será graficar n círculos concentricos de radio n*k en nuestro mapa, con las siguientes líneas de comando:

r1 <- k*1

r2 <- k*2

r3 <- k*3

n  <- seq(0,359,0.05)

a  <- seq(1,360,0.05)

theta <- (360*n)/a

     x1 <- x0 + r1*cos(theta)

     y1 <- y0 + r1*sin(theta)

     x2 <- x0 + r2*cos(theta)

     y2 <- y0 + r2*sin(theta)

     x3 <- x0 + r3*cos(theta)

     y3 <- y0 + r3*sin(theta)

circle <- data.frame(c(x1,x2,x3),c(y1,y2,y3))

N <-length(x1)                        # Tamaño del vector

color <- c()                          # Inicializo vector

color[1:N]           <- 'gold'

color[(N+1):(2*N)]   <- 'darkorange'

color[(2*N+1):(3*N)] <- 'red'

{PlotOnStaticMap(PuebMap, lat = circle[,2], 

                lon = circle[,1], cex=0.5, pch=19,

                destfile = "interm.png", 

                col= color ,add=TRUE)};

Quedará algo como esto:

Podemos definir un nuevo nivel de zoom:

Por último y con fines de análisis geoespacial, para definir en que radio se encuentra cada una de las Unidades Económicas (UE), empleamos las siguientes líneas:

D <- dim(DF)[1]; D
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 1 o no
DF$PotR1 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r1^2){
      DF[j,'PotR1'] <- 'dentro'
      } else DF[j,'PotR1'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 2 o no
DF$PotR2 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r2^2){
      DF[j,'PotR2'] <- 'dentro'
      } else DF[j,'PotR2'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 3 o no
DF$PotR3 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r3^2){
      DF[j,'PotR3'] <- 'dentro'
      } else DF[j,'PotR3'] <- 'fuera'
}

Finalmente, en este ejemplo, se analiza a las variables cuantitativas en función de su pertenencia o no al área definida por cada radio, tal como el número de posibles empleados en cada radio:

sumER_1 <- 0
for (j in 1:D){
if (DF[j,'PotR1']=='dentro'){
    sumER_1 <- DF[j,'pers']+ sumER_1
    }
}
sumER_1                       # Potencial Radio 1

sumER_2 <- 0
for (j in 1:D){
if (DF[j,'PotR2']=='dentro'){
    sumER_2 <- DF[j,'pers']+ sumER_2
    }
}
sumER_2 <- sumER_2 - sumER_1             
sumER_2                       # Potencial Radio 2

sumER_3 <- 0
for (j in 1:D){
if (DF[j,'PotR3']=='dentro'){
    sumER_3 <- DF[j,'pers']+ sumER_3
    }
}
sumER_3 <- sumER_3 - sumER_2 - sumER_1   
sumER_3                       # Potencial Radio 3


sumER_x <- sum(DF['pers']) - sumER_1 - sumER_2 - sumER_3
sumER_x                       # Potencial fuera de cualquier radio

Información que finalmente representamos en una tabla:

> potencial <- data.frame(rbind(c(sumER_1, sumER_2, sumER_3, sumER_x)

+                               )

+                         )

> colnames(potencial) <- c(' Radio 1', ' Radio 2', ' Radio 3', ' Radio 4')

> rownames(potencial) <- 'Empleados'

> round(potencial)

           Radio 1  Radio 2  Radio 3  Radio 4

Empleados     1539      372      398      128

Lo que nos indica que existen aproximadamente 1539 empleados en el radio 1, que enfatizemos es de 44 Km a partir de la coordenada central que hemos fijado.

Visualización de interacciones de Twitter con R (1 de 2)

En esta entrada muestro la visualización de interacciones de Twitter, con ayuda de R y tomando sólo como ejemplo el hashtag "#PRI", una vez que han pasado las elecciones.

Los paquetes necesarios son "graphTweets", "networkD3" y "igraph". Se emplea la información de 1,000 tweets del 5 de junio de 2016. Nuestra red lucirá así:

 Y si analizamos la formación de comunidades, nuestra red se verá de la siguiente forma:

Existen otros paquetes de R con los cuales también podemos graficar las interacciones como "rgl" o "tcltk".

Plot de la función de producción Cobb-Douglas en R

El propósito de esta entrada, es mostrarles principalmente una, de las seguramente muchas formas que existen en R, para generar un gráfico en tres dimensiones de la función de producción Cobb-Douglas (linealizada en este caso).

Cabe destacar que no se requiere de paquetería especial alguna, pero que el código necesario para producir el gráfico es ligeramente extenso y es posible que nos perdamos en algún momento. Mi recomendación es que lo copies, y que trates de "desmenuzarlo" poco a poco, con el objetivo de comprenderlo y posteriormente realizar las modificaciones que se ajusten a tus necesidades.

El primer paso es la base de datos, para los que se encuentran familiarizados con la función de producción Cobb-Douglas resultará sencillo entender la estructura que debe tener la base de datos, para los que no, recomiendo que lean un poco al respecto a fin de poder seguir todo el proceso. Los datos deben lucir de la siguiente manera:

En la imagen se muestran los primeros 6 datos de la base que utilizaremos en este ejemplo y que tiene el nombre de cd. La columna Y representa el valor de producción, K el valor del factor Capital, y finalmente L es el valor del factor Trabajo; todos estos valores son en realidad, el logaritmo neperiano del verdadero valor monetario ($), efectuado así para representar linealmente la función.

Lo siguiente que debemos hacer es estimar los coeficientes (elasticidades) de los factores de producción con el siguiente código:

cd.lm <-     lm(formula= I(Y) ~ I(K) + I(L),

             data  = cd, na.action=na.exclude );summary( cd.lm )

Call:

lm(formula = I(Y) ~ I(K) + I(L), data = cd, na.action = na.exclude)

Residuals:

    Min      1Q  Median      3Q     Max 

-0.2812 -0.0938 -0.0030  0.0904  0.3340 

Coefficients:

            Estimate Std. Error t value             Pr(>|t|)    

(Intercept)   2.7235     0.2644   10.30 < 0.0000000000000002 ***

I(K)          0.4566     0.0632    7.23       0.000000000067 ***

I(L)          0.8845     0.0577   15.32 < 0.0000000000000002 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.131 on 111 degrees of freedom

Multiple R-squared:  0.817,     Adjusted R-squared:  0.814 

F-statistic:  248 on 2 and 111 DF,  p-value: <0.0000000000000002

A partir de este punto, es que comenzamos a plotear la función de producción, si sigues paso a paso el siguiente código, podrás identificar que parte del gráfico es la que se esta desarrollando.

getOption("device") ()                    ## Abre ambiente gráfico
interval <- 1
x <-   seq(from = floor(min(cd$K)/interval)*interval,
           to   = ceiling(max(cd$K)/interval)*interval,
           by   = interval)

y <-   seq(from = floor(min(cd$L)/interval)*interval,
           to   = ceiling(max(cd$L)/interval)*interval,
           by   = interval)

f <-function(x,y){
   coef(cd.lm)[1] + x*coef(cd.lm)[2] + y*coef(cd.lm)[3]
   }

z <- outer(x,y,f)

zlim <- c(floor(min(z)), ceiling(max(z)))
#Alternativa zlim <- c(floor(min(z)), 1)

cd$Y.hat <- f(x = cd$K, y = cd$L)

## Comienza el plot 3D con argumentos específicos

cd.persp <- persp(
          x        = x,
          y        = y,
          z        = z,
          theta    = 310,
          phi      = 0,
          col      = 'cornflowerblue',
          ltheta   = -135,
          shade    = 0.01,
          ticktype = 'detailed',
          main     = 'Función de producción',
          xlab     = 'Capital',
          ylab     = 'Trabajo',
          zlab     = 'Valor producción',
          zlim     = zlim,
          scale    = TRUE,
          border   = NA,
          nticks   = 4 ,
          font.lab = 2)

# Se plotean los valores reales

cd.below <- subset(cd, cd$Y <= cd$Y.hat)

cd.below.trans3d <- trans3d(
       x    = cd.below$K,
       y    = cd.below$L,
       z    = cd.below$Y,
       pmat = cd.persp )

points( cd.below.trans3d, col = "red", pch = 16)

cd.below.hat.trans3d <- trans3d(
      x    = cd.below$K,
      y    = cd.below$L,
      z    = cd.below$Y.hat,
      pmat = cd.persp )

segments(
   x0  = cd.below.trans3d$x,
   y0  = cd.below.trans3d$y,
   x1  = cd.below.hat.trans3d$x,
   y1  = cd.below.hat.trans3d$y,
   lty = "solid",
   col = "red")

cd.above <- subset (cd, cd$Y >= cd$Y.hat)

cd.above.trans3d <- trans3d(
       x    = cd.above$K,
       y    = cd.above$L,
       z    = cd.above$Y,
       pmat = cd.persp )

points (cd.above.trans3d, col = "green", pch = 16)

cd.above.hat.trans3d <- trans3d(
        x    = cd.above$K,
        y    = cd.above$L,
        z    = cd.above$Y.hat,
        pmat = cd.persp )

segments(
     x0  = cd.above.trans3d$x,
     y0  = cd.above.trans3d$y,
     x1  = cd.above.hat.trans3d$x,
     y1  = cd.above.hat.trans3d$y,
     lty = "solid",
     col = "green")

# Agregamos la leyenda

legend("bottom",
       c("Ajustado","Real (superior)","Real (inferior)"),
       horiz=TRUE,
       pch=c(18,19,19),
       col=c("cornflowerblue","green","red"),
       inset=0.95,
       pt.cex=c(2,1.1,1.1),
       bty="n")

El resultado será (dependiendo de tus datos), un gráfico como el que se presenta al inicio. Sin embargo, tenemos otra vía que requiere muchas menos líneas de código; aunque cabe mencionar que existen algunas diferencias en cuanto a estética se refiere. Para elaborarlo requerimos de un paquete denominado scatterplot3d. Una vez instalada esta paquetería, procedemos con el siguiente código:

library(scatterplot3d)

cd.lm <-     lm(formula= I(Y) ~ I(K) + I(L),

             data  = cd, na.action=na.exclude );summary( cd.lm )

s3d <- scatterplot3d(cd$K,cd$L,cd$Y,
       pch=16,highlight.3d=TRUE,
       type="h",main="Función de producción")
s3d$plane3d(cd.lm)


Lo que produce un gráfico como el siguiente:

Claro está que existen más argumentos para esta función (scatterplot3d), lo que podría producir un gráfico de mejor calidad, aunque corresponde al usuario indagar para lograr el objetivo que busca.

Plotear un corazón 3D en R

En esta ocasión, después de no publicar en un largo tiempo, me doy a la tarea escribir algo con menor carga académica, pero que quizás en determinado momento, ya sea por curiosidad o por alguna necesidad tengas el interés en replicar e inclusive mejorar.

Plotear un corazón (de San Valentín) en dos o tres dimensiones en algún Software no es algo novedoso, sin duda que no lo es, existen numerosos códigos en Matlab, Mathematica, C, etc., regados por Internet. Pero para los que utilizamos mayormente R, aquí esta la alternativa.

El único paquete que necesitamos es misc3d, y plotear la siguiente función:

En donde, tanto x, y como z, son vectores que van desde -3 a 3, con intervalos k, una vez explicado esto, lo único que debemos hacer escribir las siguientes lineas de código:

# install.packages('misc3d')
library(misc3d)

k <- 0.030
x <- seq(-3, 3, by = k);x
y <- seq(-3, 3, by = k);y
z <- seq(-3, 3, by = k);z

f  <- function(x ,y, z)(-(x^2)*(z^3)-(9/80)*(y^2)*(z^3)) + ((x^2)+(9/4)*(y^2)+(z^2)-1)^3
contour3d(f,0,x,y,z, color='red', smooth=TRUE, fill= TRUE)

Existen algunos argumentos de la función que bien valdría la pena explorar, todo acorde a tus necesidades. También es valioso comentar que el intervalo k, entre más pequeño sea, mejor apariencia tendrá nuestro gráfico, pero el tiempo que tomará R para plotearlo será mayor.