Acerca de… NA’s, valores perdidos.

Hola

Aun recuerdo como me ilusionó saber que si dentro del paréntesis de la media pones na.rm=T el R pasa de los NA.

> v<-c(1,2,3,NA,4)
> mean(v)
[1] NA
> mean(v,na.rm=T)
[1] 2.5

Por que esto es muy cómodo, de hecho lo puedes meter en tu propia función para la media:


# datos de ejemplo
> coches<-cars
# metamos un NA
> coches[10,1]<-NA
# la media
> apply(coches,MARGIN=2,FUN=mean)
speed  dist
NA 42.98
# nuestra función
> mimedia<-function(x) mean(x, na.rm=TRUE)
> apply(cars,MARGIN=2,FUN=mimedia)
speed  dist
15.40 42.98

Y meter la función mimedia donde te de la gana.
Así podia deshacerme de los puñeteros NA que no sabía muy bien de donde venian ni para que servian.

Claro, con el tiempo uno aprende que los NA dependiendo de donde vengan pueden significar cosas muy distintas y que a veces podemos ignorarlos.
Pero a veces conviene convertirlos en un valor interesante.

coches[is.na(coches)]<-"valorinteresante"

O borrarlos:

coches<-cars
coches[10,1]<-NA
coches.sin.na<-coches[!is.na(coches$speed),]

Esto viene a cuento de esta entrada donde explica como son los NA según de donde vengan. Sabiendo eso podemos trabajar de una forma u otra. Ya que a veces la distribución de los NA puede estar sesgada y entonces no vale borrarlos o ignorarlos.

jaume

Acerca de… Jornadas de ususarios de R

Hace ya más de una semana que terminaron y creo que es interesante hablar de ellas aquí.

Ya está diponible el libro de resumenes. Está muy bien por que junto al resumen te da la opción de ir a la ponencia.

Hay varias sobre ecología:

Algunas aplicaciones de R en ecología: de la docencia elemental a la investigación avanzada.

Marcelino de la Cruz Rot. Universidad Politécnica de Madrid.

Presentación

Detección de patrones espaciales de biodiversidad de árboles y mamíferos en la Península Ibérica.

Jennifer Morales Barbero, Rafael Francisco García Vázquez y Dolores Ferrer Castán. Área de Ecología, Universidad de Salamanca.

Presentación

Correcciones taxonómicas de grandes bases de datos de forma automatizada.

Luis Cayuela. Universidad Rey Juan Carlos.

Presentación

Además hay algunos sobre supervivencia que se podrían aplicar a algunos estudios ecológicos.

jaume

Acerca de… RStudio

Hola

Dicen que rectificar es de sabios. Ya voy por el tercer intermediario (GUI) para R.

Todo empezó con el JGR, necesitaba hacer el uso de R un poco más fácil en linux (vedlo aquí y aquí)

La verdad es que no iba demasiado bien y era un lío trabajar de dos formas diferentes en linux y Windows. Traté de buscar uno que funcionara en ambas plataformas y apareció el RKWard. Va bastante bien pero en windows da algunos problemas a la hora de generar los gráficos en formatos vectoriales. Y el asunto de las actualizaciones tampoco es trivial.

Pero el otro día vi a un colega usando RStudio. Yo no dejaba el RKward por que te marca la sintaxis con colorines y te sugiere a medida que escribes. Pero el RStudio hace eso y mucho más, además funciona en Windows y Linux. Y la cuestión es que va suave, da una sensación de maquina bien engrasada. El tema de las actulizaciones tampoco es automático, pero no va mal, ya os contaré.

Pues eso probad el RStudio.

jaume

Acerca de… rkward.

Hola,

editado 27/10/2011: Mirad esta entrada.

Estoy empezando a utilizar una nueva interfaz gráfica para R en linux. Se llama RKward y va bastante mejor que JGR.

A pesar de ser para KDE va muy bien el Ubuntu, se instala directamente con el Añadir y Quitar de Ubuntu y no necesita un montón de librerias de KDE para funcionar. Además en su web en forma de wiki esta todo bastante bien explicado.

Puedes ver gráficamente lo que tienes en el Workspace,  tiene la opción de autocompletar con el tabulador, puedes buscar en la ayuda de forma gráfica, etc.

Una vista del RKward funcionando en Ubuntu.

Ya iré colgando cositas.

Saludos.

Acerca de… ¿Y si no sabemos qué especie es?

Hola,

Los ecólogos, al no ser taxónomos expertos, muchas veces nos encontramos que al trabajar en el campo no sabemos exactamente con que especies estamos tratando, faltan las flores, los frutos, no es la época adecuada, etc.

Texto alternativo

Las plantas no siempre están así, preparadas para ser determinadas.

Para reflejar esa incertidumbre en las publicaciones viene bien esta nomenclatura sobre la que siempre tuve dudas:

Nombredegenero sp. Se refiere a una especie determinada dentro de ese género, pero no sabemos cual.

Nombredegenero spp. Se refiere a unas cuantas especies de ese genero, pero no todas, si no solo pondríamos Nombredegenero.

cf. (confer) y aff. (affinis), por lo que he podido averiguar, se refieren a una identificación no segura, Nombredegenero cf./aff epíteto se refiere a que no estamos seguros de que sea de esa especie. La diferencia entre estas dos ideas tiene sentido en un contexto taxonómico: cf.  (comparar con…) quiere decir que la planta que intentamos determinar se debería comparar con pliegos, tipos, etc del taxón indicado. Mientras que aff. (emparentado con…) nos indica que la planta que tenemos entre manos esta emparentada con la que indicamos.

Un ejemplo.

Medicago cf. Medicago minima: Si encontramos un pliego de herbario con esta indicación nos está indicando que el/la que determinó la planta pensaba que era M. minima, y esperaba compararla con otros pliegos.

Medicago aff. minima: Un pliego de herbario con este nombre indica que la planta que contiene el pliego se parece a M. minima pero que el/la que la determinó no afirma que fuera M. minima.

Como veis la diferencia no es muy clara.

Por cierto estas abreviaturas no se escriben en cursiva.

Como otras veces lo busqué para mi, pero lo dejo aquí por si a alguien le sirve.

De hecho en ecología, el problema de no saber con qué especie estas tratando, te puede arruinar un diseño experimental. En el blog de Luís Cayuela hay una entrada sobre como tratar con esta incertidumbre estadísticamente (con R por supuesto).

Espero comentarios, criticas, enmiendas… y bibliografía, no he encontrado un sitio donde hable de esto claramente.

Saludos.

P.D.

Más sobre nomenclatura botánica aquí o aquí.

Acerca de… glm con datos trucados.

Siempre que he intentado entender la salida de un lm o un glm en el R me he encontrado con el mismo problema, como no conozco bien los datos, o son muy complicados, me queda la duda de si estaré interpretando bien lo que me dice el R.

Lo que hice el otro día, después de discutir con un compañero de trabajo que interpretaba las salidas de una forma diferente que yo, fue crear unos datos con los que sabía lo que iba a salir. Así por fin entendí del todo la salida del glm y me di cuenta de que interpretaba mal las salidas.

A continuación os pego lo que hice, primero como crear unos datos trucados, después una exploración con algunos gráficos y por fin los análisis.

Espero que os sea útil.

Creemos los datos trucados.

Creamos los datos de base aleatorios

origen<-rnorm(900,mean=0,sd=10)

es un vector de 900 elementos distribuidos como una normal de media 0 y desviación típica 10

Creamos los factores:

un vector de 900 elementos con 300 para cada nivel del factor

factorA<-c(rep(“controlA”,300)
,rep(“a1″,300)
,rep(“a2″,300)
)

y un vector de 900 elementos que divide cada nivel del factor anterior en 3 grupos de 100

factorB<-c(rep(
c(rep(“controlB”,100)
,rep(“b1″,100)
,rep(“b2″,100)
)
,3
)
)

Ahora los hacemos factores y arreglamos los niveles de base.

factorA<-as.factor(factorA)
factorA<-relevel(factorA,”controlA”)
#El control es el nivel base… así después en el análisis queda todo mejor
factorB<-as.factor(factorB)
factorB<-relevel(factorB,”controlB”)        #id

Hasta aquí tenemos la estructura típica de unos datos para analizar. Ahora hay que hacer el efecto de los factores sobre los datos.
Creamos los efectos
Modificación producida por cada nivel del factor… enseguida tendrá sentido.

controlA<-0
a1<-20
a2<-30
controlB<-0
b1<-6
b2<-12

Vector de los efectos, este vector modificará los datos origen.

efectosA<-c(rep(controlA,300),rep(a1,300),rep(a2,300))
efectosB<-c(rep(c(rep(controlB,100),rep(b1,100),rep(b2,100)),3))

Si os fijáis los valores están en la misma posición en que están los niveles de los factores.

Ahora creamos los datos definitivos, creamos los datos modificando el vector aleatorio con los efectos de los factores

datos<-origen + efectosA + efectosB

Si os fijáis lo que estamos haciendo aquí es coger una distribución normal y modificarla en función de unos factores que es lo que en realidad suponemos que está pasando cuando vamos a analizar unos datos haciendo un ANOVA.

si hubiera interacción, cuando se juntan el a1 y b2 el efecto se multiplicaría.

datos.int<-replace(datos,factorA==”a1″&factorB==”b2″,datos[factorA=="a1"&factorB=="b2"]*2)

Replace mola, lo que pone aarriba le dice al R: “en el vector datos, cuando el vector factorA es igual a a1 o el vector factorB es igual a b2, pon lo que haya pero multiplicado por 2″.

Y ya tenemos los datos trucados.
Con estos datos se puede jugar modificando los valores de ControlA, ControlB, a1, a2, b1 y b2.

Exploremos los datos trucados.

DATOS SIN INTERACCIÓN:

plot(datos~factorA)

plot(datos~factorB)

La diferencia entre los niveles de los factores es, aproximadamente, la que hemos dado al crear los datos.

Veamos las medias:

tapply(datos,list(factorA,factorB),mean)
controlB            b1                 b2
controlA   -0.9258813      6.894822     11.98925
a1           20.3019293      27.829253   31.74711
a2           29.1412938      36.630104   41.87341

Y un gráfico de interacción:

interaction.plot(factorA,factorB,datos,fun=mean)

Como veis es bastante soso, al pasar de un nivel a otro de A se incrementa su valor y lo mismo ocurre con B. Pero si probamos con los datos con interacción sale distinto.

DATOS CON INTERACCIÓN

plot(datos.int~factorA)
plot(datos.int~factorB)
interaction.plot(factorA,factorB,datos.int,fun=mean)

Fijaos en las dos primeras gráficas; Ha aumentado la variabilidad en los grupos en los que hemos hecho la interacción. Y en el gráfico de interacción ya está muy claro como la línea correspondiente a b2 se comporta de forma rara.

Análisis

Probemos con los datos sin interacción, por suerte los datos son normales, je je.
Primero hacemos (los puristas dirán ajustamos) el modelo.
modelo.glm<-glm(datos~factorA+factorB+factorA:factorB)
summary(modelo.glm)

Call:
glm(formula = datos ~ factorA + factorB + factorA:factorB)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-33.4490   -6.7897   -0.3912    6.3939   34.0443

Coefficients:
Estimate       Std. Error     t value     Pr(>|t|)
(Intercept)               -0.3451     1.0159         -0.340     0.734
factorAa1                 19.9125     1.4367        13.860     < 2e-16 ***
factorAa2                 28.9171    1.4367         20.128     < 2e-16 ***
factorBb1                  8.0833     1.4367          5.626     2.46e-08 ***
factorBb2                 12.7213     1.4367         8.855      < 2e-16 ***
factorAa1:factorBb1   -1.2277     2.0317        -0.604      0.546
factorAa2:factorBb1   -0.8083     2.0317        -0.398      0.691
factorAa1:factorBb2    0.1649     2.0317         0.081      0.935
factorAa2:factorBb2    1.1961     2.0317         0.589      0.556

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

(Dispersion parameter for gaussian family taken to be 103.1990)

Null deviance: 249869  on 899  degrees of freedom
Residual deviance:  91950  on 891  degrees of freedom
AIC: 6738

Number of Fisher Scoring iterations: 2

Después, hacemos el ANOVA.
anova(modelo.glm, test=”F”)

Analysis of Deviance Table

Model: gaussian, link: identity

Response: datos

Terms added sequentially (first to last)

Df    Deviance   Resid.Df   Resid. Dev        F Pr(>F)
NULL                                     899        249869
factorA             2   131624      897        118245           637.7182 <2e-16 ***
factorB             2     26171      895        92075            126.7971 <2e-16 ***
factorA:factorB  4        124      891        91950             0.3012 0.8772

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

Como se ve en el modelo ninguna de las interacciones entre los niveles de ambos factores es significativa, después el ANOVA nos lo enseña de otra forma.

Probemos con los datos con interacción.

Primero el modelo:
modelo.glm.int<-glm(datos.int~factorA+factorB+factorA:factorB)
summary(modelo.glm.int)

Call:
glm(formula = datos.int ~ factorA + factorB + factorA:factorB)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-46.2969   -6.9360   -0.3912    6.7193   44.3720

Coefficients:
Estimate      Std. Error    t value   Pr(>|t|)
(Intercept)               -0.3451       1.1532         -0.299   0.765
factorAa1                 19.9125       1.6309        12.210   < 2e-16 ***
factorAa2                 28.9171       1.6309        17.731   < 2e-16 ***
factorBb1                  8.0833       1.6309        4.956     8.60e-07 ***
factorBb2                 12.7213       1.6309       7.800     1.72e-14 ***
factorAa1:factorBb1   -1.2277       2.3064       -0.532    0.595
factorAa2:factorBb1   -0.8083       2.3064       -0.350    0.726
factorAa1:factorBb2  32.6186       2.3064      14.143    < 2e-16 ***
factorAa2:factorBb2    1.1961       2.3064        0.519    0.604

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

(Dispersion parameter for gaussian family taken to be 132.9891)

Null deviance: 432750  on 899  degrees of freedom
Residual deviance: 118493  on 891  degrees of freedom
AIC: 6966.3

Number of Fisher Scoring iterations: 2

Después el ANOVA:

anova(modelo.glm.int, test=”F”)

Analysis of Deviance Table

Model: gaussian, link: identity

Response: datos.int

Terms added sequentially (first to last)

Df    Deviance   Resid. Df    Resid. Dev       F           Pr(>F)
NULL                                    899          432750
factorA              2      176817     897         255933            664.779   < 2.2e-16 ***
factorB              2      90566       895         165367            340.502   < 2.2e-16 ***
factorA:factorB   4      46874       891         118493            88.116     < 2.2e-16 ***

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

Como se puede ver, ahora el resumen del modelo – summary() -  nos da como significativa la interacción entre a1 y b2, justo los datos que son diferentes entre datos y datos.int. Además en el ANOVA también sale significativa la interacción.

También está aquí.

Ahora habría que ver qué pasa si modificamos para 2 niveles de dos factores, pero eso para otro día.

Saludos.

Acerca de… R en el New York Times.

Hola.

He encontrado otro blog de wordpress en el que escriben sobre R. El blog se llama Modelos estadisticos 1, al parecer es el blog de un curso sobre modelos estadísticos. Pinchando aquí podéis leer una entrada en la que hacen referencia a un artículo en el New York Times sobre R.

El autor de dicho artículo,Ashlee Vance, además tiene una entrada sobre R en un blog del NYTimes. Según dice, ya hay una empresa que trabaja con una versión comercial de R, algo parecidoa RedHat en linux.

Aunque no es la única hay otra llamada Vhayu, mirad aquí para un comentario en el blog de la empresa.

Saludos.

Acerca de… trabajar con dataframes en R

Hola.

Los dataframes son una forma muy útil de organizar los datos en R. Pero siempre me encuentro con dificultades para trabajar con ellos, sobretodo para hacer cosas triviales, por ejemplo ¿Como se crea o se borra una columna?

Usemos el dataframe dades del post anterior.

Para ver una columna de un dataframe basta con usar $, hay que poner

> nombredf$nombrecolumna
> dades$datos1

teniendo en cuenta esto podemos trabajar con estos datos como si fueran un vector…

> mean(dades$datos1)
[1] 51.877

…nos da lo mismo que nos daba el summary de otro post anterior.

Pero ¿Y si ahora queremos añadir otra columna al dataframe?

Primero tenemos que asegurarnos de que la longitud del dataframe sea la misma que la de la nueva columna.

> length(dades$datos1)
[1] 40

le pido la longitud sólo de una de las columnas, por que al ser un dataframe todas tienen que medir lo mismo.

Creamos la una nueva columna de 40 elementos

> nueva.col<-c(seq(1:40))

y la pegamos al dataframe

> dades$factor3<-nueva.col

si es un factor se lo decimos

> dades$factor3<-as.factor(dades$factor3)

y ya está

> summary(dades)
factor1   factor2  datos1           datos2          datos3          factor3
a:20      1:10     Min.   : 0.32    Min.   : 12.17  Min.   : 1.00   1: 1
b:20      2:10     1st Qu.:25.35    1st Qu.:263.51  1st Qu.:10.75   2: 1
          3:10     Median :54.05    Median :476.13  Median :30.00   3: 1
          4:10     Mean   :51.88    Mean   :501.87  Mean   :30.00   4: 1
                   3rd Qu.:76.39    3rd Qu.:725.57  3rd Qu.:49.25   5: 1
                   Max.   :99.96    Max.   :996.24  Max.   :59.00   6:1
                                                                    (Other):34

Para eliminarla es un poco raro. Sabemos que para ver un elemento de un vector usamos los corchetes

> nueva.col[3]
[1] 3

nos muestra el elemento nº3 del vector nueva.col

para ver una columna, fila o celda de un dataframe los corchetes funcionan parecido [fila,columna], para ver el contenido de la segunda columna:

> dades[,2]
[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4
[39] 4 4
Levels: 1 2 3 4

Si queremos borrar una columna tenemos que usar el signo “-”, por ejemplo, para borrar la columna factor3 (que es la 6ª) del dataframe:

> dades[,-6]

Esto nunca me ha parecido muy elegante, pero no he encontrado otra forma de hacerlo. Si alguien sabe una forma mejor que lo diga, o mejor, que lo escriba.

Saludos.

Acerca de… más de JGR

Hola.

editado 27/10/2011: Mirad esta entrada.

Más cosas sobre JGR:

Se puede hacer un ejecutable en el escritorio para ejecutar JGR directmente sin ejecutar el R.

Nos lo sugiere en la consola cuando ponemos JGR().

Recordemos, lo primero cargar la libreria:

> library(JGR)
Loading required package: rJava
Loading required package: JavaGD
Loading required package: iplots

Please use the corresponding JGR launcher to start JGR.
Run JGR() for details. You can also use JGR(update=TRUE) to update JGR.

Aquí ya nos lo dice, pero no nos dice donde está el launcher, pero si ponemos:
> JGR()

nos dice:

Starting JGR …
(You can use /usr/local/lib/R/site-library/JGR/scripts/run to start JGR directly)

Supongo que según la estructura de directorios de cada uno será diferente.

Lo que hay que hacer es crear un lanzador en el escritorio (botón derecho sobre el escritorio – crear lanzador) y poner ese path donde pone comando: desde ese momento nos saldrá el icono en el escritorio y al pincharlo se ejecutará el JGR directamente.

Otra cosa, he tenido problemas para que funcione la ayuda, para que funcione correctamente la ayuda hay que instalar además el paquete r-base-html. No es un paquete de R como el base o el graphics que se instala desde R. Es un paquete de Ubuntu que se instala con:

sudo apt-get install r-base-html

o desde Synaptic.

Espero que os haya sido util.

Sobre como instalar el JGR mirad aquí. He editado aquella entrada para que esté esto último también.

Por cierto, también hay un JGR para windos.

Saludos.

Acerca de… cargar datos en R

Hola.

Una explicación sobre como cargar datos en R a partir de hojas de cálculo. Hay más formas, pero yo suelo usar esta.

He creado en calc una hoja de cálculo con 5 columnas llamada dades.ods, tiene dos columnas llamadas factor1 y factor2 y tres columnas con los datos.

pantalazo dades.ods

Leamos esto con el R:

Se puede leer directamente de .ods o .xls usando read.dbf, pero Nunca sabemos de donde van a venir los datos así que lo explico para leerlos desde un archivo de texto .csv que se puede crear desde cualquier plataforma, o casi.

Pasamos dades.ods a .csv mediante Guardar como (también se puede hacer con MSExcel)

Antes de leerlo con R miradlo con un editor de texto para ver que la separación por campos está bien hecha, el texto está entre comillas, los decimales con puntos, etc. Si no hacemos bien esto último y los decimales son comas, el R pensará que ahí hay un cambio de campo y nos haremos un lío.

Una vez comprobado esto leemos el archivo dades.csv desde R:

dades->read.table("dades.csv",sep=",",header=T)

le decimos que bajo el nombre dades nos ponga el contenido de dades.csv, que los campos estan separados por comas (sep=”,”) y que nuestro archivo tenía cabeceras (header=T). Al crear el .csv podemos especificar esto.

Como veis no he especificado ruta, supongo que dades.csv está en el directorio de trabajo.

Para ver el contenido:

> summary(dades)
 factor1    factor2         datos1          datos2           datos3
 a:20    Min.   :1.00   Min.   : 0.32   Min.   : 12.17   Min.   : 1.00
 b:20    1st Qu.:1.75   1st Qu.:25.35   1st Qu.:263.51   1st Qu.:10.75
         Median :2.50   Median :54.05   Median :476.13   Median :30.00
         Mean   :2.50   Mean   :51.88   Mean   :501.87   Mean   :30.00
         3rd Qu.:3.25   3rd Qu.:76.39   3rd Qu.:725.57   3rd Qu.:49.25
         Max.   :4.00   Max.   :99.96   Max.   :996.24   Max.   :59.00

o

> str(dades)
'data.frame':   40 obs. of  5 variables:
 $ factor1: Factor w/ 2 levels "a","b": 1 1 1 1 1 1 1 1 1 1 ...
 $ factor2: int  1 1 1 1 1 2 2 2 2 2 ...
 $ datos1 : num  62.04 0.82 73.25 60.96 52.31 ...
 $ datos2 : num  504 996 725 259 338 ...
 $ datos3 : int  1 2 3 4 5 6 7 8 9 10 ...

Nos ha creado un dataframe llamado dades.

Ahora hay que decirle qué columnas son factores y cuáles son datos, ya que los que eran texto (factor1) sí que los ha reconocido como factores, pero los que eran números (factor2) nos los ha reconocido como número y en el summary nos ha hecho la estadística.

dades$factor2<-as.factor(dades$factor2)

Le decimos que cree un factor nuevo llamado dades$factor2 que será el mismo dades$factor2 pero convertido en factor.

Si hacemos ahora el summary:

> summary(dades)
 factor1 factor2     datos1          datos2           datos3
 a:20    1:10    Min.   : 0.32   Min.   : 12.17   Min.   : 1.00
 b:20    2:10    1st Qu.:25.35   1st Qu.:263.51   1st Qu.:10.75
         3:10    Median :54.05   Median :476.13   Median :30.00
         4:10    Mean   :51.88   Mean   :501.87   Mean   :30.00
                 3rd Qu.:76.39   3rd Qu.:725.57   3rd Qu.:49.25
                 Max.   :99.96   Max.   :996.24   Max.   :59.00

tenemos tres variables datos 1-3 y dos factores: factor1 con dos niveles a y b y factor2 con 4 niveles 1, 2, 3, 4.

Saludos.