Introducción a R

• Descargar el instalador de https://www.python.org/downloads/windows, por ejemplo, Latest Python 3 Release

• Una vez instalado, ejecutar en el terminal cmd.exe (quizá haya que pulsar intro alguna vez):

  py -m pip install --user --upgrade pip
  py -m pip install --user numpy scipy pandas statsmodels matplotlib sklearn 
  
  

• Otra posibilidad es usar conda
  • sistema de paquetería basado en Python
  • ventajas:
    • multiplataforma
    • facilita distribución de programas (p.ej. qiime2)
    • no requiere privilegios de administrador
  • inconveniente: duplica muchas bibliotecas instaladas por la paquetería del sistema operativo

apt install python3-pandas python3-dill python3-sklearn # módulos necesarios para tratar datos
apt install python-numpy-doc python-scipy-doc python-pandas-doc python-sklearn-doc # documentación

Conviene instalar algún editor, por ejemplo,

apt install idle-python3.8
apt install spyder
apt install emacs
apt install jupyter jupyter-notebook python3-ipykernel

• Descargar el último R desde https://cran.r-project.org/bin/windows/base
• Recomendaciones si se va a usar Rcmdr
  • Tras ejecutar el fichero .exe descargado, pulsar Siguiente hasta el diálogo Opciones de configuración. Ahí elegir Sí y pulsar Siguiente.
  • Elegir en Modo de display la opción SDI - ventanas separadas.
  • Mantener las opciones por omisión en el resto de diálogos, hasta Instalando.
  • Al pinchar en el icono se abre el programa R Gui.
• Interfases opcionales:

  • Para facilitar el uso básico y el aprendizaje: el paquete Rcmdr (commander); instalarlo desde el menú Paquetes del RGui o escribiendo en la R Console la siguiente orden:

    install.packages ("Rcmdr")
    
    

    • Ejecutar en la consola de R la orden

      library(Rcmdr)
      
      

    • Cuando pregunte si instalar paquetes adicionales, pulsar Sí.
    • Es posible que pregunte en algún momento si instalar paquetes en un directorio alternativo. Aceptar el que da por omisión.
  • El entorno más popular: RStudio
  • Emacs Speaks Statistics desde Emacs para Windows

• Menú Datos, submenú Conjunto de datos en paquetes.
• Eligiendo alguno de los paquetes (datasets, carData…) se listan los conjuntos de datos disponibles.
• Los datos del paquete datasets (por ejemplo, mtcars) están disponibles al arrancar R pero, al cargar Rcmdr, se ocultan (para no abarrotar la lista de datos disponibles en la barra de elementos activos) y es necesario cargarlos explícitamente.
• Ejercicio
  • Cargar mtcars del paquete datasets.
  • Cargar Chile del paquete carData.

• Menú Datos, opción Cargar conjunto de datos.
• Habitualmente, dichos ficheros tienen extensión .rda ó .RData.
• Pueden contener
  • un solo conjunto de datos (cuyo nombre puede coincidir o no con el nombre del fichero)
  • múltiples conjuntos de datos, además de variables, funciones, etc. (entorno de trabajo de R)
• Ejercicio
  • Cargar mtcoches.rda y comprobar la lista de datos disponibles
  • Cargar iris.RData y comprobar la lista de datos disponibles

• Menú Datos, submenú Importar datos.
• De especial interés es la importación de datos de texto (CSV)
  • un renglón = un registro
  • campos delimitados por un carácter
    • csv = coma «,»
    • csv2 = punto y coma «;»
    • delim = tabulador
    • table = espacio en blanco (espacios a tabuladores; uno o varios)
  • codificación: utf-8 (unicode) o latin-1 (iso-8859-1)
• Ejercicio: intenta importar los siguientes conjuntos de datos
  • vacas.csv
  • vacas.txt
  • vacas.blancos1
  • vacas.blancos2

Para encontrar ayuda sobre un comando de R, por ejemplo rep, puedes usar

?rep

Para algunos comandos con nombres especiales, hay que usar comillas; por ejemplo,

?"["
?"if"

Si sólo conocemos parte del nombre del comando (p.ej. norm), usamos

apropos("norm")

Si queremos encontrar información relacionada con un concepto, podemos usar:

??variance
help.search("variance")

• https://cran.r-project.org/doc/contrib/R-intro-1.1.0-espanol.1.pdf
• «Spanish» en https://cran.r-project.org/other-docs.html
• tutoriales sobre técnicas variadas: quick R
• paquetes agrupados por campos: https://cran.r-project.org/web/views/
• gráficos: R graph gallery

help("regression")
help("modules regression")

• Dive into Python 3
• Python Data Science Handbook

• condicionales

  Generemos un vector de temperaturas:

  t <- runif (15, 10, 40)
  
  

  Para combinar comparaciones vectorialmente, usamos ifelse y un operador de circuito largo («&» = «y»; «|» = «ó»):

  ifelse (t>20 & t<30, "agradable", "desagradable")
  
  

  Para combinar comparaciones de forma que esperamos un único booleano, usamos operadores de circuito corto («&&» = «y»; «||» = «ó»):

  n=10; if (is.numeric(n) && length(n)==1 && n>0) sqrt(n) else "no puedo"
  
  

• bloques de instrucciones mediante llaves «{}»; se usan en contextos donde tiene que ir más de una instrucción:

  if (n > 5) {print("demasiado grande"); n<-5} else {print("demasiado pequeño"); n<-n+1}
  
  

  El valor devuelto por un bloque es el valor devuelto por su última instrucción.

• bucles

  for (i in 1:10) cat ("El cuadrado de", i, "es", i^2, "\n")
  
  

  Detrás de «in» puede ir cualquier vector o lista; por ejemplo

  for (i in letters) {cat (i); cat (toupper (i)); cat ("-")}; cat ("\n")
  for (i in iris) print (summary (i))
  
  

• bucles con acumulador

  Se define una variable (acumulador) antes del bucle para que recoja los resultados:

  suma <- 0; for (i in 1:10) suma <- suma+i; suma                       # equivale a sum(1:10) -> suma
  producto <- 1; for (i in 1:10) producto <- producto*i; producto       # prod(1:10); factorial(10) -> producto
  bector <- c(); for (i in 1:10) bector <- c(bector,i^2); bector        # (1:10) ^ 2 -> bector
  lista <- list(); for (i in 1:10) lista <- c(lista,list(1:i)); lista   # sapply (1:10, function (i) 1:i) -> lista
  

• condicionales

  Para combinar comparaciones de forma que esperamos un único booleano, usamos operadores de circuito corto («and» = «y»; «or» = «ó»):

  import math
  n=10
  if ((isinstance(n,int) or isinstance(n,float)) and n>0):
     math.sqrt(n) 
  else:
     "no puedo"
  
  

  Generemos un vector de temperaturas:

  import scipy.stats as st
  t = st.uniform.rvs (10, 40, 15)
  
  

  Para combinar comparaciones vectorialmente, usamos if así:

  ["agradable" if ti>20 and ti<30 else "desagradable" for ti in t]
  
  

• nótense los bloques de instrucciones tras dos puntos y con sangría:

  if n > 5:
      print("demasiado grande")
      n = 5
  else:
      print("demasiado pequeño")
      n++
  
  

• bucles

  for i in range(10):
    print ("El cuadrado de", i, "es", i**2)
  
  mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/mtcars.csv")
  for v in mtcars:
    mtcars[v].describe()
  
  

• bucles con acumulador

  Se define una variable (acumulador) antes del bucle para que recoja los resultados:

  • numéricos

    suma = 0
    for i in range(1,11):
      suma += i
    
    

    producto = 1
    for i in range(1,11):
      producto *= i
    
    

  • recolección

    vector = []
    for i in range(1,11):
       vector.append (i**2)
    
    

    lista = []
    for i in range(10):
       lista.append (list (range (i)))
    
    

• con nombre

  f <- function(x) x^2
  f(15)
  
  

• anónimas (útiles para definiciones al vuelo, como se ve en la sección siguiente)

  (function(x) x^2) (15)
  
  

• pueden tener varios argumentos/parámetros y varias instrucciones en el cuerpo

  var.ponderada <- function (x, pesos)     # varianza ponderada
  {  mp <- weighted.mean (x, pesos)        # media ponderada
     weighted.mean ((x - mp) ^ 2, pesos) }
  
  

• cada argumento puede tener un valor por omisión

  var.ponderada <- function (x, pesos = rep(1,length(x)))
  {  mp <- weighted.mean (x, pesos)
     weighted.mean ((x - mp) ^ 2, pesos) }
  
  

• pueden pasar sus argumentos a otras funciones mediante «…»

  var.ponderada <- function (x, pesos, ...)
  {  mp <- weighted.mean (x, pesos, ...)
     weighted.mean ((x - mp) ^ 2, pesos, ...) }
  
  

  Por ejemplo:

  muestra <- c (1, 20, NA, 5); pesos <- c (1,2,3,4)
  weighted.mean (muestra, pesos); 
  weighted.mean (muestra, pesos, na.rm=TRUE)
  var.ponderada (muestra, pesos)                  # produce NA
  var.ponderada (muestra, pesos, na.rm=TRUE)      # pasa «na.rm=TRUE» a weighted.mean
  
  

• recursivas: una función puede llamarse a sí misma

  mifactorial <- function(n) if (n==0) 1 else n*mifactorial(n-1)
  
  

  otro ejemplo:

  nombres <- c("Señor don Manuel Montenegro", "Asun Lubiano", "Doña Antonia Salas",
               "Don Norberto Corral", "Carlos Carleos")
  títulos <- c("don", "doña", "señor", "señora")
  
  ## dado un vector de nombres, obtener el primer nombre de pila de cada uno
  nombre.de.pila <- function(cadena)
      nombre.de.pila_vector (strsplit (cadena, " ") [[1]]) # convierte "cadena" en vector de palabras
  
  nombre.de.pila_vector <- function(v)    # v = vector de cadenas
  {
      palabra1 <- v[1]
      if (tolower(palabra1) %in% títulos) # tolower = convierte a minúsculas
        nombre.de.pila_vector(v[-1])      # llamada recursiva
      else palabra1
  }
  
  for (nombre in nombres) print (nombre.de.pila (nombre))
  

• con nombre

  def f (x):
    return x**2
  f(15)
  
  

• anónimas (útiles para definiciones al vuelo, como se ve en la sección siguiente)

  (lambda x: x**2) (15)
  
  

  sólo pueden contener una expresión;

• pueden tener varios argumentos

  def media_ponderada (x, pesos):
     suma = 0
     for i in range (len (x)):
        suma += x[i] * pesos[i]
     return suma / sum(pesos)
  
  

• cada argumento puede tener un valor por omisión

  def media_ponderada (x, pesos = None):
     suma = 0
     for i in range (len (x)):
        suma += x[i] * (pesos[i] if pesos else 1)
     return suma / (sum(pesos) if pesos else len(x))
  
  

• pueden pasar sus argumentos a otras funciones mediante un argumento con «*»

  def var_ponderada (x, *otros):
    mp        = media_ponderada (x, *otros)
    cuadrados = [(xi-mp)**2 for xi in x]
    return media_ponderada (cuadrados, *otros)
  
  

  Por ejemplo:

  muestra = [1, 20, 100, 5]; pesos = [1,2,3,4]
  media_ponderada (muestra)
  media_ponderada (muestra, pesos)
  var_ponderada (muestra)                  # pasa *otros como vacío
  var_ponderada (muestra, pesos)           # pasa *otros como «pesos»
  
  

• sobre secuencias (vectores o listas)

  • sapply (con s de simple) se usa para generar un vector (o una matriz, si el resultado de la función aplicada es vectorial)

    sapply (1:10, sqrt)                       # aplicar una función con nombre; igual que sqrt(1:10)
    sapply (1:10, function(x) x**2)           # aplicar una función anónima; igual que (function(x) x**2) (1:10)
    sapply (1:10, function(x) c(x^2,x^3))     # el resultado es una matriz; distinto de (function(x) c(x^2,x^3)) (1:10)
    
    

  • sapply conserva etiquetas si se aplica a un vector de cadenas

    sapply(3:5, sqrt)                         # sin etiquetas
    sapply(c("uno","dos","tres"), nchar)      # con etiquetas
    structure (sapply(3:5, sqrt), names=3:5)  # forzar etiquetas
    
    

  • lapply (con l de lista) devuelve una lista

    lapply (1:10, function(i) list (i, list(letters[i],month.name[i])))
    
    

    aunque si se quieren conservar etiquetas mejor usar sapply o (si hace falta) sapply(..., simplify=FALSE)

    lapply (c ("uno", "dos", "tres"), nchar)
    sapply (c ("uno", "dos", "tres"), nchar)
    sapply (c ("uno", "dos", "tres"), nchar, simplify=FALSE)
    
    

• apply se usa para aplicar una función a filas o columnas de una matriz…

  m <- matrix (1:12, 3)
  apply (m, 1, sum)                 # aplicar por filas
  apply (m, 2, sum)                 # aplicar por columnas
  
  

  …o a un dataframe; en este caso, téngase en cuentra que se trasforma en matriz (homogénea) antes de aplicarse la función, por lo que (por ejemplo) si el dataframe contiene un factor todos los datos aparecen como cadenas:

  ## queremos trasformar todos los 4 de ChickWeight a NA
  ChickWeight
  cuatroNA <- function (x) {x[x==4] <- NA; x}
  cuatroNA (1:10)                                       # funciona
  cuatroNA (as.character (1:10))                        # funciona
  cuatroNA (factor (1:10))                              # funciona
  summary (apply (ChickWeight, 2, cuatroNA))            # apply no funciona (trasforma todo a factores o cadenas)
  summary (data.frame (lapply (ChickWeight, cuatroNA))) # lapply sí funciona
  

  es decir, cuando queremos aplicar una función a cada variable de un dataframe,

  • es mejor considerar el dataframe como lista y usar lapply
  • que considerar el dataframe como matriz y usar apply(, 2,)

  otro ejemplo usando la función identity; identity(x) devuelve x

  mtcars
  apply(mtcars, 1, identity) # queda traspuesta
  apply(mtcars, 2, identity)
  identical(mtcars, data.frame (apply (mtcars, 2, identity))) # TRUE
  m    <- mtcars
  m$am <- factor(m$am)
  apply(m, 1, identity) # genera matrices de cadenas
  apply(m, 2, identity) # idem
  identical(m, data.frame (apply (m, 2, identity))) # FALSE
  

• para aplicar por grupos:

  prod <- 1:20; sexo <- rep(c("M","H"),each=10)
  tapply (prod, sexo, mean)           # por niveles de "sexo"
  by (prod, sexo, mean)               # ídem
  
  

• la más general (mapeo)

  mapply (rep, 1:4, 4:1)              # aplica funciones con varios argumentos
  
  

• un comando muy habitual en simulaciones:

  replicate (10000, {m <- rnorm(5); mean(m)})
  
  

  Funciona parecido a

  sapply (1:10000, function (i) {m <- rnorm(5); mean(m)})
  
  

  donde la variable i es sintácticamente necesaria pero su valor no se usa.

• trasformación por descripción (mapeo)

  from math import sqrt
  [sqrt(i) for i in range(10)]
  list (map (sqrt, range(10)))               # lo mismo
  list (map (lambda x: x**2, range(10)))     # con función anónima
  
  

  otra versión de la media ponderada anteriormente definida:

  def media_ponderada (x, pesos):
     return  sum (map (lambda xi,pi: xi*pi, x, pesos)) / sum(pesos)
  
  

  otra versión con argumento opcional:

  def media_ponderada (x, pesos = None):
     return  sum (map (lambda xi,pi: xi*pi, x, [1]*len(x) if pesos is None else pesos)) / (sum(pesos) if pesos else len(x))
  
  

• aplicar una función a filas o columnas de un dataframe

  from sklearn import datasets
  iris = datasets.load_iris()        # datos de ejemplo
  import pandas as pd
  I = pd.DataFrame(iris.data, columns=iris.feature_names)
  I.apply(sum)                       # aplicar «sum» por columnas; equivale a I.sum()
  I.apply(sum,axis=0)                # ídem
  I.apply(sum,axis=1)                # aplicar «sum» por filas; equivale a I.sum(axis=1)
  
  

• para aplicar por grupos:

  I["especie"] = pd.Series(pd.Categorical.from_codes(codes=iris.target,categories=iris.target_names))
  I.groupby("especie").apply(lambda x: x.iloc[:,0:4].median()) # 0:4 para evitar "especie"
  
  

(Lo mismo en Rcmdr, Python.)

• Guardar todo el entorno de trabajo (funciones, dataframes…)

  save.image ()                                 # lo guarda en ./.RData (el fichero cargado al arrancar R)
  save.image (file = "~/trabajo/mi-sesión.rda") # en fichero arbitrario
  
  

• Guardar uno o varios objetos en formato nativo de R

  save (mi.dataframe, mi.función, file = "~/trabajo/mis-objetos.rda")
  
  

• Guardar un dataframe como una tabla de texto (CSV)

  write.table (mtcars, "/tmp/mtcars1.csv") # sep=" " dec="." sin campo de rownames en cabecera
  write.csv   (mtcars, "/tmp/mtcars2.csv") # sep="," dec="."
  write.csv2  (mtcars, "/tmp/mtcars3.csv") # sep=";" dec=","
  
  

(Lo mismo en R, Rcmdr.)

• Guardar entorno de trabajo

  import dill
  dill.dump_session ("/tmp/miTrabajo.pkl") # recuperar con dill.load_session
  
  

• Guardar un objeto en formato nativo de Python

  misDatos.to_pickle ("/tmp/misDatos.pkl")
  
  

• Guardar un dataframe como una tabla de texto (CSV)

  misDatos.to_csv ("/tmp/misDatos.csv")
  
  

(Lo mismo en Rcmdr, Python.)

D    <- mtcars                # dataframe de ejemplo
D$am <- factor (mtcars$am)    # una cualitativa; también: D$am <- factor(D$am, labels=c("auto","manual"))
summary (D)                   # resumen del dataframe
T <- table (D$cyl)            # frecuencias absolutas
prop.table (T)                # frecuencias relativas
summary (D$hp)                # media y cinco cuartiles
mean (D$hp)                   # media
var (D$hp)                    # varianza
sd (D$hp)                     # desvío típico
median (D$hp)                 # mediana
quantile (D$hp, 0.95)         # percentil 95
quantile (D$hp, 0:4/4)        # cinco cuartiles; similar a fivenum(D$hp)
IQR (D$hp)                    # recorrido intercuartílico

En muchas funciones se puede incluir la opción na.rm=TRUE para omitir los ausentes:

mean (c (1:10, NA))
mean (c (1:10, NA), na.rm=TRUE)

Gráficas:

barplot (table (D$cyl))   # diagrama de barras
pie (table (D$am))        # diagrama de sectores
hist (D$hp)               # histograma
boxplot (hp ~ am, D)      # diagrama de cajas
plot (mpg ~ hp, D)        # diagrama de dispersión (nube de puntos)
pdf ("fichero.pdf")       # guarda gráficas subsiguientes en 1 fichero PDF (vectorial)
for (i in 1:ncol(mtcars)) {
  if (is.numeric (mtcars[[i]]))
    hist (mtcars[[i]])
  else 
    barplot (table (mtcars[[i]])) }
dev.off ()                # cierra el PDF
png ("fichero%d.png")     # guarda gráficas subsiguientes en ficheros PNG (pixelados)
for (i in 1:3)
  boxplot (mtcars[[i]] ~ mtcars$am)
dev.off ()                # cierra el PNG

Ejercicio:

• Escribir un programa en R que realice gráficas adecuadas para las variables del conjunto de datos «Chile» del paquete «car».

  data (Chile, package="car")
  Chile$education <- factor (Chile$education, levels=c("P","S","PS"))
  
  

(Lo mismo en R, Rcmdr.)

Con los datos clásicos de flores iris:

import pandas as pd
i = pd.read_csv ("datos/iris.csv")
i.describe()                           # numérico
i["Species"].describe()                # moda
iSn = i['Species'].value_counts()      # recuentos
iSn / iSn.sum()                        # frecuencias relativas
i["Petal.Length"].mean()               # media
i["Petal.Length"].var()                # varianza
i["Petal.Length"].std()                # desvío típico
i["Petal.Length"].median()             # mediana
i["Petal.Length"].quantile(0.5)        # ídem
i["Petal.Length"].quantile([.05,.95])  # percentiles 5 y 95
import numpy as np                     # para «nan»
pd.Series([1,2,np.nan,4]).mean()       # excluye los ausentes

Gráficas:

import matplotlib.pyplot as plt
i["Petal.Length"].plot(kind="hist")                # histograma
plt.show()                                         # mostrar en pantalla
i["Petal.Length"].value_counts().plot(kind="bar")  # diagrama de barras
plt.savefig("barras.png")                          # guardar en fichero
i["Petal.Length"].value_counts().plot(kind="pie")  # diagrama de sectores
i["Petal.Length"].plot(kind="box")                 # diagrama de caja
i.boxplot("Petal.Length","Species")                # cajas por grupo

(Lo mismo en Python, Rcmdr.)

• Crear nueva variable de un dataframe:

  mtcars$consumo <- (3.8 / mtcars$mpg * 1.6) * 100
  
  

• Recodificar:

  mtcars$am <- factor (mtcars$am, labels=c("auto","manual"))
  mtcars$consumo3 <- cut (mtcars$consumo, c(-Inf, 10, 20, 30, Inf), c("muy poco", "poco", "mucho", "demasiado"))
  
  

• Reordenar niveles de factor:

  data (Chile, package="car")
  Chile$educación <- factor (Chile$education, levels = c("P","S","PS")) # P=primaria, S=secundaria, PS=postsecundaria
  
  

• Filtrar

  edades.mujeres <- Chile$age [Chile$sex == "F"]
  mtcars4v <- mtcars [mtcars$cyl==4 & mtcars$vs==0, ]
  ?==                           # ayuda sobre < > == != 
  ?&                            # ayuda sobre ! & |
  
  

(Lo mismo en R, Rcmdr.)

import pandas as pd
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/Chile.csv")

• Crear nueva variable en un dataframe:

  Chile.euro = 1.2 * Chile.income
  
  

• Recodificar:

  Chile.sexo = Chile.sex.astype("category")
  Chile.sexo.values.categories = ["mujer", "varón"]
  Chile.clase = pd.cut (Chile.income, [-np.inf,30000,100000,np.inf], labels=["baja","media","alta"])
  
  

• Reordenar niveles de categórica:

  Chile["educación"] = Chile["education"].astype("category",categories=["P","S","PS"])
  
  

• Filtrar:

  edades_mujeres = Chile.age[Chile.sex=="F"]
  varones_jovenes = Chile[(Chile.sex=="M") & (Chile.age<35)]
  
  

Más en la documentación de Pandas.

(Lo mismo en Python, Rcmdr.)

• letra inicial
  • r para generar aleatorios
  • p para la función de distribución \(F_X\,(t) = \Pr\,[X\leq t]\)
  • d para la función de densidad de variables continuas o la función de probabilidad (o de cuantía, o de masa) de variables discretas \(f_X\,(t) = \Pr\,[X=t]\)
  • q para la función cuantil
• abreviatura del nombre de la familia de la distribución
  • unif uniforme
  • norm gausiana
  • exp exponencial
  • binom binomial
  • pois Poisson
  • t Student
  • chisq \(\chi^2\) (ji-cuadrado)

runif (10, 0, 1)                   # genera 10 aleatorios con distribución uniforme entre 0 y 1
pnorm (150, 170, 15)               # probabilidad de que una gausiana N(170;15) sea menor que 150
plot (x <- seq(0,5,.1), dexp(x,1)) # curva de densidad de una exponencial con media 1
qpois (0.5, 1.7)                   # mediana de una distribución de Poisson con media 1,7

(Lo mismo en R, Rcmdr.)

• métodos
  • rvs para generar aleatorios
  • cdf función de distribución (cumulative density function)
  • pdf función de densidad (probability density function)
  • pmf función masa de probabilidad (probability mass function)
  • mean, var, std, median - media, varianza, desvío típico, mediana
  • ppf función cuantil (percent point function)
• distribuciones
  • uniform uniforme
  • norm gausiana
  • expon exponencial
  • binom binomial
  • pois Poisson
  • t Student
  • chi2 ji-cuadrado (\(\chi^2\))

import scipy.stats as st
st.uniform.rvs (2, 0.1, 30)     # 30 valores uniformes entre 2 y 2+0,1
st.norm.cdf (200, 170, 10)      # función de distribución en 200 de una N(170;10)
import matplotlib.pyplot as plt # función de densidad (de la docstring de expon)
fig, ax = plt.subplots(1, 1)
x = np.linspace(expon.ppf(0.01),
                expon.ppf(0.99), 100)
ax.plot(x, expon.pdf(x),
        'r-', lw=5, alpha=0.6, label='exponencial')
## también se pueden fijar previamente los parámetros:
v = st.binom (10, .3)           # binomial n=10, p=0,3
v.mean(), v.var(), v.std()      # media, varianza, desvío típico
v = st.poisson (100)            # puasón con media 100
v.ppf([.05,.95])                # percentiles 5 y 95

(Lo mismo en Python, Rcmdr.)

?cars                            # velocidad y distancia de frenada
library (car)                    # para datos Chile
## contraste de gausianidad
shapiro.test (cars$speed)        # comprobar mediante stem o hist
shapiro.test (Chile$statusquo)
shapiro.test (Chile$age)
## información del objeto "contraste"
shapiro.test (cars$speed) -> contraste
names (contraste)                # "statistic" "p.value"   "method"    "data.name"
contraste$p.value                # nivel crítico o p-valor
## simulación
n <- 100                                                        # tamaño muestral
m <- rnorm(n); stem(m); shapiro.test(m)$p.value                 # relación entre forma y p-valor
mean (replicate (10000, shapiro.test(rnorm(n))$p.value) < 0.05) # proporción de rechazos
## contraste sobre la media
t.test (cars$speed, mu=15)                        # bilateral
t.test (cars$speed, mu=15, alternative="less")    # H0: mu=15; H1: mu<15
t.test (cars$speed, mu=15, alternative="greater") # H0: mu=15; H1: mu>15
## contraste sobre proporción
summary (Chile$sex)
?prop.test
prop.test (sum(Chile$sex=="F"), nrow(Chile), 0.5) # H0: p = 1:2; H1: p =/= 1:2
## contraste sobre la mediana
wilcox.test (cars$speed, mu=15)                        # bilateral
wilcox.test (cars$speed, mu=15, alternative="less")    # H0: mu=15; H1: mu<15
wilcox.test (cars$speed, mu=15, alternative="greater") # H0: mu=15; H1: mu>15
M  <- rexp (25)                                        # exponencial con media=1
Me <- qexp (.5)                                        # mediana teórica = 0,693.147.2
wilcox.test (M, mu=Me)                                 # poco recomendable porque exige simetría
mean (replicate (100000, {M <- rexp(25); wilcox.test(M,mu=Me)$p.value<.05})) # 11% en vez de 5%
prop.test (sum(M>Me), 25, .5)                          # prueba de signos: poco potente, pero cumple
mean (replicate (100000, {M <- rexp(25); prop.test(sum(M>Me),25,.5)$p.value<.05})) # <5%
mean (replicate (100000, {M <- rexp(25); t.test(M,mu=1)$p.value<.05})) # 7,5% (pero con media, no mediana)
## contraste sobre dos proporciones
## ¿es significativamente mayor la proporción de mujeres que votan Y?
## p = proporción que votan Y; H0: pF=pM; H1: pF>pM
mujeres <- Chile [Chile$sex=="F" & !is.na(Chile$vote),]
varones <- Chile [Chile$sex=="M" & !is.na(Chile$vote),]
prop.test (c (sum (mujeres$vote=="Y"), sum (varones$vote=="Y")),
           c (nrow (mujeres),          nrow (varones)),
           alternative = "greater")
## contraste de igualdad de varianzas (auxiliar para el de igualdad de medias)
? iris                  # comparar varianza de Petal.Length en «versicolor» y «virginica»
with (iris, plot (Petal.Length, Petal.Width, col=Species))
boxplot (Petal.Length ~ Species, iris)
with (iris, by (Petal.Length, Petal.Species, shapiro.test))
PLversicolor <- iris [iris$Species=="versicolor", "Petal.Length"]
PLvirginica  <- iris [iris$Species=="virginica",  "Petal.Length"]
var.test (PLversicolor, PLvirginica)      # contraste de igualdad de varianzas
var.test      (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # lo mismo
bartlett.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # otro similar
## contraste de igualdad de medias (muestras independientes)
t.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica")) # estimando dos varianzas distintas
t.test (Petal.Length ~ Species, iris, subset = Species %in% c("versicolor","virginica"),
        var.equal = TRUE)                                                                # casi no cambia el p-valor
## contraste de igualdad de medias (muestras pareadas)
? ChickWeight # ¿Hay incremento significativo de peso a partir del día 20?
boxplot (weight ~ Time, ChickWeight)
cw <- ChickWeight [ChickWeight $ Time %in% c (20, 21), ]
t.test (weight ~ Time, cw, alternative="less") # como si fueran independientes
cw20       <- cw [cw$Time == 20,]; rownames(cw20) <- cw20$Chick
cw21       <- cw [cw$Time == 21,]; rownames(cw21) <- cw21$Chick
ambos      <- intersect (rownames(cw20), rownames(cw21))
incremento <- cw21[ambos,"weight"] - cw20[ambos,"weight"]
shapiro.test (incremento) # innecesario por el tamaño de muestra (45)
t.test (cw20[ambos,"weight"], cw21[ambos,"weight"], paired=TRUE, alternative="less")
t.test (incremento, mu=0, alternative="greater") # equivale al anterior
## contraste sobre la mediana de la diferencia (muestras independientes)
## comparar Petal.Width en «setosa» y «virginica» dentro del conjunto de datos «iris»
?iris
shapiro.test (iris [iris$Species == "setosa",    "Petal.Width"])
shapiro.test (iris [iris$Species == "versicolor","Petal.Width"])
t.test      (Petal.Width ~ Species, iris, subset = Species %in% c("setosa","versicolor")) # suficiente tamaño muestral
wilcox.test (Petal.Width ~ Species, iris, subset = Species %in% c("setosa","versicolor")) # resultado parecido
## contraste sobre la mediana de la diferencia (muestras pareadas)
?sleep            # ¿hay diferencias significativas en sueño extra según somnífero?
boxplot (extra ~ group, sleep)
extra1 <- sleep$extra [sleep$group == '1']
extra2 <- sleep$extra [sleep$group == '2']
dif    <- extra1 - extra2
shapiro.test (dif)
wilcox.test (extra1, extra2, paired=TRUE)       # asegurarse de que el orden se corresponde en ambos vectores
wilcox.test (sleep ~ group, sleep, paired=TRUE) # otra forma de expresarlo
wilcox.test (dif)                               # es equivalente a un contraste de Wilcoxon de 1 muestra
t.test (extra1, extra2, paired=TRUE)            # no muy distinto

• Otro ejemplo:

  • Comparar las dos secuencias «x» y «y» de seudoaleatorios de cierto generador. Cada valor de «y» es el consecutivo (dentro de la secuencia seudoaleatoria) del correspondiente en «x».

    ?randu                                             
    dif <- randu$x - randu$y                           # comparar con el siguiente seudoaleatorio
    shapiro.test (dif)                                 # irrelevante por el tamaño muestral, en principio...
    wilcox.test (randu$x, randu$y, paired=TRUE)
    t.test (randu$x, randu$y, paired=TRUE)             # ...pero p-valor cerca del 5%
    ## comparar con seudoaleatorios decentes
    xy <- matrix(runif(400*5),400,byrow=TRUE)[,1:2]; x <- xy[,1]; y <- xy[,2]
    wilcox.test (x, y, paired=TRUE)
    

(Lo mismo en R, Rcmdr.)

import numpy  as np                # para nan, isnan
import pandas as pd                # para dataframes
from scipy import stats            # para contrastes
import matplotlib.pyplot as plt    # para histograma
cars  = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/cars.csv") 
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/Chile.csv")
## contraste de gausianidad
stats.shapiro (cars.speed)
stats.shapiro (cars.speed) [1]                     # p-valor
stats.shapiro (Chile.age [~ np.isnan (Chile.age)]) # hay NA
## simulación
n = 100                                                     # tamaño muestral
## relación entre forma y p-valor
m = stats.norm.rvs(0,1,n); m.plot(kind="hist"); plt.show(); stats.shapiro(m)[1]
(pd.Series([stats.shapiro(stats.norm.rvs(0,1,n))[1] for i in range(10)]) < 0.05) . mean() # proporción de rechazos
## contraste sobre la media
pvalor2 = stats.ttest_1samp (cars.speed, 15) . pvalue                       # bilateral
## ya que mediamuestral = 15,4 > 15 entonces
pvalor1 = pvalor2 / 2                                                       # unilateral H0: mu=15; H1: mu>15
1 - pvalor1                                                                 # unilateral H0: mu=15; H1: mu<15
## contraste sobre proporción H0: p = 1:2; H1: p =/= 1:2
Chile.sex.value_counts()
n = Chile.sex.count()
p = Chile.sex.value_counts(True)['F']
from math import sqrt
2 * (1 - stats.norm.cdf (abs ((p - 0.5) / sqrt (0.25/n))))  # p-valor bilateral
stats.chisquare ([1379,1321]) . pvalue                      # otra posibilidad: usar contraste ji-2 equivalente
## contraste sobre la mediana
M  = stats.expon.rvs (0, 1, 25)                        # exponencial con media=1
Me = stats.expon.ppf (.5)                              # mediana teórica = 0,693.147.2
stats.wilcoxon (M - Me)                                # poco recomendable porque exige simetría
nM = sum (M-Me > 0); stats.chisquare ([nM, len(M)-nM]) # prueba de signos: poco potente
## contraste sobre dos proporciones
## ¿es significativamente mayor la proporción de mujeres que votan Y?
## p = proporción que votan Y; H0: pF=pM; H1: pF>pM
tabla = Chile.groupby(['sex','vote']).size().unstack()                       # tabla de contingencia
stats.chi2_contingency(tabla)[1]                                             # P-valor de contraste ji-2 de homogeneidad
mujeres = Chile [(Chile.sex=="F") & ~pd.isnull(Chile.vote)]
varones = Chile [(Chile.sex=="M") & ~pd.isnull(Chile.vote)]
pvalor2 = stats.chi2_contingency([[sum(varones.vote=="Y"),sum(varones.vote!="Y")],
                                  [sum(mujeres.vote=="Y"),sum(mujeres.vote!="Y")]]) [1]
pF = sum(mujeres.vote=="Y")/len(mujeres) # 37%
pM = sum(varones.vote=="Y")/len(varones) # 32%
pvalor1 = pvalor2 / 2                    # porque pF > pM
## contraste de igualdad de varianzas (auxiliar para el de igualdad de medias)
## ¿hay diferencias significativas en sueño extra entre los dos grupos?
sleep = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/sleep.csv")
stats.shapiro (sleep.extra [sleep.group == 1])
stats.shapiro (sleep.extra [sleep.group == 2])
stats.bartlett (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2])
## contraste de igualdad de medias (muestras independientes)
stats.ttest_ind (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2]) # suponiendo varianzas iguales
stats.ttest_ind (sleep.extra [sleep.group == 1], sleep.extra [sleep.group == 2], equal_var=False)
## contraste de igualdad de medias (muestras pareadas)
## ¿hay incremento significativo de peso a partir del día 20?
CW = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/ChickWeight.csv")
CW.boxplot ("weight", "Time")
plt.show()
CW.weight[CW.Time==20].mean()                                              # menor
CW.weight[CW.Time==21].mean()                                              # mayor
stats.ttest_ind (CW.weight[CW.Time==20], CW.weight[CW.Time==21]).pvalue/2  # como si fueran independientes
cw = CW[CW.Time.isin([20,21])]
cw20 = cw [cw.Time == 20]; cw20.index = cw20.Chick
cw21 = cw [cw.Time == 21]; cw21.index = cw21.Chick
ambos = cw20.index.intersection(cw21.index)
stats.shapiro (cw20.loc[ambos,"weight"] - cw21.loc[ambos,"weight"]) # innecesario por el tamaño de muestra (45)
stats.ttest_rel (cw20.loc[ambos,"weight"], cw21.loc[ambos,"weight"]).pvalue / 2
## contraste sobre la mediana de la diferencia (muestras independientes)
## comparar Petal.Width en «setosa» y «virginica» dentro del conjunto de datos «iris»
iris = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/iris.csv")
stats.shapiro (iris.loc[iris.Species == "setosa",    "Petal.Width"])
stats.shapiro (iris.loc[iris.Species == "virginica", "Petal.Width"])
stats.ttest_ind (iris.loc[iris.Species=="setosa","Petal.Width"], 
                 iris.loc[iris.Species=="virginica","Petal.Width"], equal_var=False)
stats.mannwhitneyu (iris.loc[iris.Species=="setosa","Petal.Width"], 
                    iris.loc[iris.Species=="virginica","Petal.Width"], alternative='two-sided') # parecido
## contraste sobre la mediana de la diferencia (muestras pareadas)
randu = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/randu.csv") # comparar «x» con «y»
dif = randu.x - randu.y
stats.shapiro (dif)                   # irrelevante por el tamaño muestral
stats.wilcoxon (dif)                 
stats.wilcoxon  (randu.x, randu.y)    # otra forma de hacer lo mismo
stats.ttest_rel (randu.x, randu.y)    # parecido
stats.ttest_1samp (dif, 0)            # lo mismo

(Lo mismo en Python, Rcmdr.)

library(car) # para Chile

### contraste ji-2

## ¿hay relación entre región y sexo?
table (Chile$region, Chile$sex)
prop.table(table (Chile$region, Chile$sex), 1)
chisq.test (table (Chile$region, Chile$sex))
## ¿hay relación entre educación y sexo?
## ¿hay relación entre sexo y voto?
## ¿hay relación entre educación y voto?

### contraste de correlación lineal

## ¿hay relación entre población, edad, ingresos y statusquo?
variables <- c("population", "age", "income", "statusquo") # names (which (sapply (Chile, is.numeric)))
cor    (Chile[variables], use="pair")
sapply (Chile[variables], shapiro.test)                    # no gausianas => mejor Spearman
plot (Chile$population, Chile$age)
cor.test (Chile$population, Chile$age)                     # por omisión, Pearson (exige gausianidad)
cor.test (Chile$population, Chile$age, method="spearman")
plot (Chile$population, Chile$income)
plot (jitter(Chile$population), jitter(Chile$income))      # jitter ayuda para representar discretas con pocos niveles
## etc.

(Lo mismo en R, Rcmdr.)

import numpy as np                 # para isnan
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt    # para scatter
Chile = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/Chile.csv")

### contraste ji-2

## ¿hay relación entre región y sexo?
tabla = Chile.groupby(["region","sex"]).size().unstack()
tabla = pd.crosstab (Chile.region, Chile.sex)               # otra forma
pd.crosstab (Chile.region, Chile.sex, normalize="index")    # index=filas columns=columnas
stats.chi2_contingency (tabla)
## ¿hay relación entre educación y sexo?
## ¿hay relación entre sexo y voto?
## ¿hay relación entre educación y voto?

### contraste de correlación lineal

## ¿hay relación entre población, edad, ingresos y statusquo?
Chile.corr()
variables = ["population", "age", "income", "statusquo"]
Chile[variables].apply(stats.shapiro)                                  # pero hay «nan»es
Chile[variables].apply(lambda x:stats.shapiro(x[~np.isnan(x)]))        # no gausianas => mejor Spearman
plt.scatter (Chile.population, Chile.age)
def cor_test (x, y):
    xsinnan = ~ np.isnan (x)
    ysinnan = ~ np.isnan (y)
    ambos   = xsinnan & ysinnan
    return stats.pearsonr (x[ambos], y[ambos])
cor_test (Chile.population, Chile.age)                                 # pearson
stats.spearmanr (Chile.population, Chile.age)                          # spearman
## etc.

(Lo mismo en Python, Rcmdr.)

cor (mtcars)                       # buscamos predictor para «mpg»
sort (abs (cor(mtcars)[,"mpg"]))   # el peso parece el mejor predictor lineal
plot (mpg ~ wt, mtcars)
modelo <- lm (mpg ~ wt, mtcars)
abline (modelo, col=2, lwd=3)      # añade recta; lwd = linewidth = anchura
summary (modelo)                   # coeficientes, R2, p-valor, estima de sigma
## p-valor del modelo = p-valor de la pendiente
predict (modelo, data.frame(wt=3)) # para un coche de 3.000 libras
coef(modelo) %*% c(1,3)            # lo mismo; 1 multiplica al intercepto

(Lo mismo en R, Rcmdr.)

import pandas as pd
mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/mtcars.csv")
mtcars.corr()
abs(mtcars.corr().mpg).sort_values() # «wt» parece el mejor predictor lineal
import matplotlib.pyplot as plt
x = mtcars.wt;              y = mtcars.mpg
x = x.values.reshape(-1,1); y = y.values.reshape(-1,1)
plt.scatter (x, y)
plt.show()
from sklearn import linear_model
regre = linear_model.LinearRegression()
regre.fit (x, y)
plt.scatter (x, y)
plt.plot (x, regre.predict (x), color='red', linewidth=3)
plt.show()
regre.predict(3)                     # predicción para 3.000 libras
regre.intercept_ + regre.coef_ * 3   # lo mismo
regre.score (x, y)                   # R2

Con otro módulo (quizá requiera «apt install python-statsmodels»):

import pandas as pd
mtcars = pd.read_csv ("http://bellman.ciencias.uniovi.es/cursoR/dat/mtcars.csv")
import statsmodels.formula.api as smf
regre = smf.ols ("mpg ~ wt", data = mtcars)
resul = regre.fit()
print (resul.summary ())
print (resul.summary2 ()) # incluye estimación de varianza residual (scale)

representación

• redondeo

  floor(pi)                               # por abajo
  ceiling(pi)                             # por arriba
  round(pi)                               # redondeo a entero
  round(pi,3)                             # redondeo a milésimas
  round(max(mtcars$hp),-2)                # redondeo a centenas
  signif(c(pi,1/3), 3)                    # cifras significativas
  formatC(5, width=10, flag="0")          # forzar anchura; rellenar con ceros
  

• coma decimal

  options()$OutDec                        # por omisión, punto
  options("OutDec")
  options(OutDec=",")
  1/3
  options(OutDec="'")
  1e-3
  

representación en gráficas

• abbreviate: recorta automáticamente de forma inambigua

  > abbreviate(rownames(mtcars))
            Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive 
               "MRX4"              "MRXW"              "D710"              "Hr4D" 
    Hornet Sportabout             Valiant          Duster 360           Merc 240D 
               "HrnS"              "Vlnt"              "D360"              "M240" 
             Merc 230            Merc 280           Merc 280C          Merc 450SE 
               "M230"             "Mr280"             "M280C"            "M450SE" 
           Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental 
            "Mr450SL"           "M450SLC"              "CdlF"              "LncC" 
    Chrysler Imperial            Fiat 128         Honda Civic      Toyota Corolla 
               "ChrI"              "F128"              "HndC"            "TytCrl" 
        Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28 
             "TytCrn"              "DdgC"              "AMCJ"              "CZ28" 
     Pontiac Firebird           Fiat X1-9       Porsche 914-2        Lotus Europa 
               "PntF"              "FX1-"              "P914"              "LtsE" 
       Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
               "FrPL"              "FrrD"              "MsrB"              "V142" 
  

• matrices multidimensionales (array)

  is.array(Titanic) # TRUE
  dim(Titanic)      # 4 2 2 2
  dimnames(Titanic) # lista
  T <- Titanic
  dimnames(T) <- list(Clase=c("1ª","2ª","3ª","Tripul."),
                      Sexo=c("Varón","Mujer"),
                      Edad=c("Niño","Adulto"),
                      Sobrevivió=c("No","Sí"))
  
  

• un array es un vector con atributo dim; otra forma de definir una matriz:

  > a <- 1:12
  > attributes(a)
  NULL
  > attr(a,"dim") <- c(3,4)
  > a
       [,1] [,2] [,3] [,4]
  [1,]    1    4    7   10
  [2,]    2    5    8   11
  [3,]    3    6    9   12
  

• también se puede usar apply con tensores

  h <- HairEyeColor
  dim(h)                                  # 4 4 2
  dimnames(h)
  apply (h, "Sex", sum)                   # pueden usarse cadenas en vez de índices
  apply (h, c("Hair","Eye"), sum)         # equivale a: apply (h, 1:2, sum)
  

• ejercicio
  • calcula el porcentaje de supervivientes según la clase en el tensor Titanic

• se pueden pegar verticalmente con rbind

  d <- data.frame (x = 1:5, y = c("uno","dos","tres","cuatro","cinco"), z = 11:15)
  rbind (d[3:5,], d[1:2,])
  
  

  pero pueden dar problemas si tienen distintas columnas:

  (d1 <- d[1:2])
  (d2 <- d[2:3])
  rbind(d1,d2)                            # Error in match.names(clabs, names(xi)) : 
  					#   names do not match previous names
  
  mi.rbind <- function (x, y) # para dataframes
  {
      nx <- names(x)
      ny <- names(y)
      for (nombre in setdiff (nx, ny))
  	y[[nombre]] <- rep(NA,nrow(y))
      for (nombre in setdiff (ny, nx))
  	x[[nombre]] <- rep(NA,nrow(x))
      rbind (x, y)
  }
  mi.rbind(d1,d2)
  ## o recurriendo a biblotecas:
  dplyr::bind_rows(d1,d2)
  

• se pueden cruzar usando merge indicando las columnas con identificadores

  D <- data.frame(id = c(1,4,5), grupo = c("A","A","B"))
  merge (d, D, by.x="x", by.y="id")
  merge (d, D, by.x="x", by.y="id", all=TRUE)
  
  

  esta última versión incluye todos los registros, rellenando con NA en caso necesario

• ejercicios
  • combina en un único dataframe los datos 2010..2018/vale002.csv de ./dat/aireUO/dat/csv

  • combina en un único dataframe

    • ped y dat
    • ped y dat0

    del fichero ./dat/datped.rda

• ejemplo
  • sea \(X\sim\mathcal B(n=10,p=0.7)\) y sea una muestra x <- rbinom (100, 10, 0.7)
  • sea \( Y_i \sim\mathcal P(\lambda=x_i) \) y sea y <- rpois (100, x)

• secuencias

  plot (x)
  plot (x, type = "p") # puntos
  plot (x, type = "l") # líneas
  plot (x, type = "b") # ambos (both)
  plot (x, type = "h") # «histograma»
  
  

• gráfico de dispersión de y sobre x

  plot (x, y)
  plot (y ~ x)               # «fórmula»
  plot (y ~ x, datos)        # si «x» y «y» en dataframe «datos»
  with (datos, plot (y ~ x)) # idem
  
  

• gráfico de una función en un intervalo

  plot (sin, -pi, pi)
  
  

• cajas

  boxplot(x)
  
  

• histograma

  hist(x)
  
  

• tallos y hojas

  stem(x)
  
  

• barras y sectores

  barplot(table(x))
  pie(table(x))
  
  

Para añadir elementos a gráficos ya creados mediante plot, etc.:

• recta \(y=a+b\cdot x\)

  abline(a,b)
  
  

  Por ejemplo:

  plot (y ~ x)
  abline (lm (y ~ x), col=2, lwd=3) # col = color; lwd = anchura de la recta
  
  

• lines (x, y) para añadir líneas
• points (x, y) para añadir puntos
• text (x, y, etiquetas) para añadir textos

Para que los gráficos vayan a un fichero pixelado. Por ejemplo:

png()
plot(y~x)
plot(x~y)
dev.off()    # para que cierre correctamente los ficheros .png

Para que los gráficos vaya a un fichero vectorial.

pdf()              # mete todos los gráficos subsiguientes en Rplots.pdf
pdf(onefile=FALSE) # crea un .pdf por cada gráfica

Hay que terminar usando dev.off() al igual que con png()

Sirve para que los gráficos se cierren correctamente. Si no se usase, podrían quedar operaciones pendiente y el gráfico podría estar sin terminar.

Para guardar un gráfico ya expuesto en una ventana.

savePlot()

• predecir Y a partir de X: \(Y = \beta_0 + \beta_1\,X + \epsilon\)
• los \(\beta\) se estiman a partir de los datos \(x_i,y_i\), \(i=1,...,n\) minimizando los residuos \(e_i\): \[\sum_{i=1}^n e_i^2 = \sum \left(y_i - \beta_0-\beta_1\,x_i \right)^2 \]
• la calidad predictiva suele expresarse mediante el coeficiente de determinación \[ R^2 = \frac {\mbox{variación de Y explicada por el modelo}} {\mbox{variación total de Y}} \]

• para construir un modelo, elegir la variable más relacionada con la respuesta

  > sort (cor (mtcars) [,1] ^ 2)
       qsec      gear      carb        am        vs      drat        hp      disp 
  0.1752963 0.2306734 0.3035184 0.3597989 0.4409477 0.4639952 0.6024373 0.7183433 
        cyl        wt       mpg 
  0.7261800 0.7528328 1.0000000 
  > summary (lm (mpg ~ wt, mtcars))
  
  Call:
  lm(formula = mpg ~ wt, data = mtcars)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -4.5432 -2.3647 -0.1252  1.4096  6.8727 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
  wt           -5.3445     0.5591  -9.559 1.29e-10 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 3.046 on 30 degrees of freedom
  Multiple R-squared:  0.7528,Adjusted R-squared:  0.7446 
  F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10
  

• predecir con intervalo de confianza:

  > predict (modelo, data.frame (wt = 5), interval="confidence", level=0.99)
         fit      lwr      upr
  1 10.56277 7.447369 13.67817
  > predict (modelo, data.frame (wt = 5), interval="confidence", level=0.95)
         fit     lwr      upr
  1 10.56277 8.24913 12.87641
  
  > abscisas <- data.frame (wt = seq (min(mtcars$wt), max(mtcars$wt), length.out=100))
  > p <- predict (modelo, abscisas, interval="prediction")
  > head (p)
         fit      lwr      upr
  1 29.19894 22.58903 35.80885
  2 28.98781 22.39104 35.58458
  3 28.77667 22.19276 35.36059
  4 28.56554 21.99420 35.13688
  5 28.35441 21.79535 34.91346
  6 28.14327 21.59621 34.69033
  > plot (mpg ~ wt, mtcars)
  > abline (modelo <- lm (mpg ~ wt, mtcars), col = 2)
  > lines (abscisas$wt, p[,"lwr"], col=3)
  > lines (abscisas$wt, p[,"upr"], col=3)
  

• diagnosis

  • ¿se ajusta la nube de puntos a la recta?: diagrama de dispersión

    plot (X, Y)
    
    

  • residuos

    plot (lm (Y ~ X))
    
    

    • deben no presentar estructura
    • el gráfico QQ (cuantil - cuantil) muestra su gausianidad

• trasformación a regresión lineal

  plot (1/mpg ~ wt, mtcars)
  lm (1/mpg ~ wt, mtcars) -> modelo
  summary (modelo)
  abline (modelo, col = 2)
  abscisas <- data.frame (wt = seq (min(mtcars$wt), max(mtcars$wt), length.out=100))
  p <- predict (modelo, abscisas, interval="prediction")
  lines (abscisas$wt, p[,"lwr"], col=3)
  lines (abscisas$wt, p[,"upr"], col=3)
  

  probar también mpg ~ log(wt)

• para expresiones más complejas puede usarse «nls»

  > x <- -(1:100)/10
  > y <- 100 + 10 * exp(x / 2) + rnorm(x)/10
  > nlmod <- nls(y ~  Const + A * exp(B * x))
  Warning message:
  In nls(y ~ Const + A * exp(B * x)) :
    No starting values specified for some parameters.
  Initializing ‘Const’, ‘A’, ‘B’ to '1.'.
  Consider specifying 'start' or using a selfStart model
  > nlmod
  Nonlinear regression model
    model: y ~ Const + A * exp(B * x)
     data: parent.frame()
    Const       A       B 
  99.9807 10.0175  0.4964 
   residual sum-of-squares: 0.8255
  
  Number of iterations to convergence: 10 
  Achieved convergence tolerance: 2.925e-07
  

• lenguaje de modelos de R: capítulo 11 de Una introducción a R

• una variable y trasformaciones suyas

  > modelo2 <- lm (mpg ~ wt + I(wt^2), mtcars)
  > summary (modelo2)
  
  Call:
  lm(formula = mpg ~ wt + I(wt^2), data = mtcars)
  
  Residuals:
     Min     1Q Median     3Q    Max 
  -3.483 -1.998 -0.773  1.462  6.238 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  49.9308     4.2113  11.856 1.21e-12 ***
  wt          -13.3803     2.5140  -5.322 1.04e-05 ***
  I(wt^2)       1.1711     0.3594   3.258  0.00286 ** 
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 2.651 on 29 degrees of freedom
  Multiple R-squared:  0.8191,Adjusted R-squared:  0.8066 
  F-statistic: 65.64 on 2 and 29 DF,  p-value: 1.715e-11
  > abscisas <- data.frame (wt = seq (min(mtcars$wt), max(mtcars$wt), length.out=100))
  > p <- predict (modelo2, abscisas, interval="prediction")
  > plot (mpg ~ wt, mtcars)
  > lines (abscisas$wt, p[,"fit"], col=2)
  > lines (abscisas$wt, p[,"lwr"], col=3)
  > lines (abscisas$wt, p[,"upr"], col=3)
  

• varias variables

  > summary (lm (mpg ~ ., mtcars))   # «.» significa «el resto de columnas»
  
  Call:
  lm(formula = mpg ~ ., data = mtcars)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -3.4506 -1.6044 -0.1196  1.2193  4.6271 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
  (Intercept) 12.30337   18.71788   0.657   0.5181  
  cyl         -0.11144    1.04502  -0.107   0.9161  
  disp         0.01334    0.01786   0.747   0.4635  
  hp          -0.02148    0.02177  -0.987   0.3350  
  drat         0.78711    1.63537   0.481   0.6353  
  wt          -3.71530    1.89441  -1.961   0.0633 .
  qsec         0.82104    0.73084   1.123   0.2739  
  vs           0.31776    2.10451   0.151   0.8814  
  am           2.52023    2.05665   1.225   0.2340  
  gear         0.65541    1.49326   0.439   0.6652  
  carb        -0.19942    0.82875  -0.241   0.8122  
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 2.65 on 21 degrees of freedom
  Multiple R-squared:  0.869,Adjusted R-squared:  0.8066 
  F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07
  
  > step (lm (mpg ~ ., mtcars)) -> mejor.modelo
  > summary (mejor.modelo)
  
  Call:
  lm(formula = mpg ~ wt + qsec + am, data = mtcars)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -3.4811 -1.5555 -0.7257  1.4110  4.6610 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   9.6178     6.9596   1.382 0.177915    
  wt           -3.9165     0.7112  -5.507 6.95e-06 ***
  qsec          1.2259     0.2887   4.247 0.000216 ***
  am            2.9358     1.4109   2.081 0.046716 *  
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 2.459 on 28 degrees of freedom
  Multiple R-squared:  0.8497,Adjusted R-squared:  0.8336 
  F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11
  

• EJERCICIO
  • Aplicar «step» al modelo «lm (mpg ~ wt + I(wt²), mtcars)» para obtener un modelo nuevo.
  • Comparar gráficamente predicciones e intervalos de ambos modelos.

Se encierran las instrucciones dentro de una llamada a system.time:

> system.time ({Sys.sleep(5); Sys.sleep(5)})
   user  system elapsed 
   0.00    0.00   10.01 

Se llama a proc.time antes y después de la ejecución, y se restan los valores:

comienzo <- proc.time ()
Sys.sleep(5)
Sys.sleep(5)
final <- proc.time ()
final - comienzo

Sobre todo a la hora de realizar bucles, podemos encontrar molesto que la ocurrencia de un error en una iteración implique la pérdida de los cálculos de las demás iteraciones. Para evitarlo se puede usar try :

sapply (c(3,2,1,0,1,2), function(x) if (x) 10*x else no.existo)
sapply (c(3,2,1,0,1,2), function(x) try (if (x) 10*x else no.existo))

Para aprovechar los resultados, puede ser necesario usar «as.numeric»:

as.numeric (sapply (c(3,2,1,0,1,2), function(x) try (if (x) 10*x else no.existo)))

o usar inherits(...,"try-error") como se explica en la ayuda de try

Filter (function (x) !inherits(x,"try-error"),
        lapply (c(3,2,1,0,1,2),
                function(x) try (if (x) 1:x else no.existo,
                                 silent = TRUE)))

Si en la definición de una función incluimos una llamada a browser(), entramos en el depurador y podemos:

• evaluar expresiones, incluyendo variables locales
• ejecutar paso a paso con la orden n (next)
• continuar con la ejecución mediante c
• interrumpir la ejecución mediante Q (quit)

He aquí un ejemplo de sesión, en la que se define una función para calcular el coefciente de variación:

> coef.var <- function (x) {media <- mean(x); dt <- sd(x); dt/media}
> muestra <- c(0, -9, 5, 0, 4)
> coef.var (muestra)
[1] Inf
> coef.var <- function (x) {browser(); media <- mean(x); dt <- sd(x); dt/media}
> coef.var (muestra)
Called from: coef.var(muestra)
Browse[1]> debug en #1: media <- mean(x)
Browse[2]> n
debug en #1: dt <- sd(x)
Browse[2]> media
[1] 0
Browse[2]> n
debug en #1: dt/media
Browse[2]> dt
[1] 5.522681
Browse[2]> c
[1] Inf
> 

Hemos explorado qué valores toman dt y media dentro de la función, de forma que comprobamos in situ que se realiza una división por cero.

Si queremos evitar avisos que prevemos, podemos usar suppressWarnings :

edad <- c("15", "62", "no desvela su edad", "48", "no contesta")
as.numeric (edad)
suppressWarnings (as.numeric (edad))

Puede interesar también la opción «silent=TRUE» de try.

Para usar ciertas funciones es aconsejable cambiar el comportamiento habitual de R con los NA (datos ausentes). Por ejemplo, podemos tener interés en incorporar a los datos originales las puntuaciones (scores) resultantes de un anásisis de componentes principales:

options ("na.action")       # na.omit por omisión

## ejemplo de ?na.omit ?na.exclude
DF <- data.frame(x = c(1, 2, 3), y = c(0, 10, NA))
na.omit(DF)
p <- princomp (~ x + y, DF) # análisis de componentes principales
p$scores                    # sólo para dos individuos

options (na.action = na.exclude)
p <- princomp (~ x + y, DF) # ACP, lo mismo que antes
p$scores                    # también para el tercero (aunque tenga NA)
DF$acp1 <- p$scores[,1]     # lo que permite pegarlo al dataframe original

La interfaz Rcmdr efectúa esa acción de forma predeterminada.

La asignación de opciones puede hacerse de varias formas:

options (na.action   = na.exclude)   # se asigna una función
options (na.action   = "na.exclude") # se asigna una cadena con el nombre de la función
options ("na.action" = na.exclude)   # el nombre de la opción puede estar entrecomillado
options ("na.action" = "na.exclude")

El resultado de cada asignación puede comprobarse con

options ("na.action")
options()$na.action       # mismo efecto

pero

options (na.action)

da error, porque los argumentos de "options" han de ser cadenas (con el nombre de las opciones) o parejas de la forma: etiqueta = valor.

grep ("[aeiou]",  letters) # las vocales
grep ("[^aeiou]", letters) # las consonantes (el circunflejo niega)
grep ("[a-dx-z]", letters) # las cuatro primeras letras o las tres últimas
grep ("^V.*a$", v)         # empieza por «V» y acaba en «a»

ficheros <- system ("ls -a ~", intern=TRUE) # todos los de la carpeta personal
grep ("^\\.", ficheros, value=TRUE)         # ficheros ocultos; empiezan por .

En este último caso,

• grep ("^.", ...) no funcionaría porque . se refiere a cualquier carácter
• grep ("^.", ..., fixed=TRUE) no funcionaría porque fixed=TRUE impediría el significado especial de ^
• grep (".", ..., fixed=TRUE) no funcionaría porque incluiría ficheros no ocultos cuyo nombre contiene un punto
• grep ("^\.", ...) no funcionaría porque \ es el «carácter de escape» de las cadenas de R, y \. no es una secuencia de escape reconocida

• grep ("^\\.", ...) funciona porque la cadena ^\\. es la representación externa en R de la cadena interna ^\.:

  > nchar ("^\\.")
  [1] 3
  > substr ("^\\.", 1,1)
  [1] "^"
  > substr ("^\\.", 2,2)
  [1] "\\"
  > substr ("^\\.", 3,3)
  [1] "."