4.7: Respuestas a los ejercicios

Respuesta a la pregunta de shapiro.test () estructura de salida. Primero, tenemos que recordar que casi todo lo que vemos en la consola R, es el resultado de print () 'ing algunas listas. Para extraer el componente de una lista, podemos llamarlo por signo de dólar y nombre, o por corchetes y número (si el componente no se nombra). Comprobemos la estructura con str ():

Figura\(\PageIndex{1}\) Filotaxis. De izquierda a derecha: hojas dispuestas por fórmulas 1/2, 1/3 y 2/5 de filotaxis.

Código\(\PageIndex{1}\) (R):

str(shapiro.test(rnorm(100)))

Bueno, lo más probable es que el valor p provenga del componente p.value, esto es fácil. Compruébalo:

Código\(\PageIndex{2}\) (R):

set.seed(1683)
shapiro.test(rnorm(100))$p.value

Esto es lo que queremos. Ahora podemos insertarlo en el cuerpo de nuestra función.

Respuesta al ejercicio de “normalidad del abedul”. Primero, necesitamos verificar los datos y entender su estructura, por ejemplo con url.show (). Entonces podemos leerlo en R, verificar sus variables y aplicar la función Normality () a todas las columnas apropiadas:

Código\(\PageIndex{3}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
Normality <- function(a) {ifelse(shapiro.test(a)$p.value < 0.05, "NOT NORMAL", "NORMAL")}
str(betula) # asmisc.r
sapply(betula[, c(2:4, 7:8)], Normality) # asmisc.r

(Obsérvese cómo solo se seleccionaron columnas no categóricas para la comprobación de normalidad. Usamos Str () porque ayuda a verificar números de variables, y muestra que dos variables, LOBES y WINGS tienen datos faltantes. No hay problema en usar str () en su lugar.)

Solo CATKIN (longitud del catkin hembra) está disponible para los métodos paramétricos aquí. Es un caso frecuente en datos biológicos.

¿Qué pasa con la comprobación gráfica para la normalidad, histograma o gráfica QQ? Sí, debería funcionar pero tenemos que repetirlo 5 veces. Sin embargo, el paquete de celosía permite hacerlo en dos pasos y encajar en una gráfica de enrejado (Figura\(\PageIndex{2}\)):

Código\(\PageIndex{4}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
betula.s <- stack(betula[, c(2:4, 7:8)])
qqmath(~ values | ind, data=betula.s, panel=function(x) {panel.qqmathline(x); panel.qqmath(x)})

(La celosía de biblioteca requiere un formato de datos largo donde todas las columnas se apilan en una sola y los datos se suministran con la columna de identificador, es por eso que utilizamos la función stack () y la interfaz

Hay muchas parcelas de enrejado. Por favor, compruebe usted mismo el histograma de enrejado:

Código\(\PageIndex{5}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
betula.s <- stack(betula[, c(2:4, 7:8)])
bwtheme <- standard.theme("pdf", color=FALSE)
histogram(~ values | ind, data=betula.s, par.settings=bwtheme)

(También hubo un ejemplo de cómo aplicar el tema en escala de grises a estas parcelas).

Figura Gráficas de enrejado QQ de\(\PageIndex{2}\) normalidad para las cinco variables de medición en el conjunto de datos de betula (las variables deben leerse de abajo hacia arriba).

Como se puede ver, SCALE.L también podría aceptarse como “aproximadamente normal”. Entre otros, LEAF.MAXW es “menos normal”.

Respuesta al ejercicio de variabilidad de caracteres de abedul. Para crear una función, es bueno comenzar desde prototipo:

Código\(\PageIndex{6}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
CV <- function(x) {}

Este prototipo no hace nada, pero en el siguiente paso puedes mejorarlo, por ejemplo, con el comando fix (CV). Después pruebe CV () con algún argumento sencillo. Si el resultado no es satisfactorio, vuelva a fijar (CV). Al final de este proceso, su función (en realidad, “envuelve” el cálculo de CV explicado anteriormente) podría verse así:

Código\(\PageIndex{7}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
CV <- function(x)
{
100*sd(x, na.rm=TRUE)/mean(x, na.rm=TRUE)
}

Entonces sapply () podría usarse para verificar la variabilidad de cada columna de medición:

Código\(\PageIndex{8}\) (R):

betula <- read.table("http://ashipunov.info/shipunov/open/betula.txt", h=TRUE)
CV <- function(x){100*sd(x, na.rm=TRUE)/mean(x, na.rm=TRUE)}
sapply(betula[, c(2:4, 7:8)], CV)

Como se puede ver, LEAF.MAXW (ubicación del ancho máximo de hoja) tiene la mayor variabilidad. En el asmisc.r, hay la función CVs () que implementa esta y otras tres mediciones de variación relativa.

Responde a la pregunta sobre los datos dact.txt. El archivo complementario dact_c.txt lo describe como un extracto aleatorio de algunas mediciones de plantas. Desde el primer capítulo, sabemos que es sólo una secuencia de números. En consecuencia, scan () sería mejor que read.table (). Primero, cargue y verifique:

Código\(\PageIndex{9}\) (R):

dact <- scan("data/dact.txt")
str(dact)

Ahora, podemos comprobar la normalidad con nuestra nueva función:

Código\(\PageIndex{10}\) (R):

Normality <- function(a) {ifelse(shapiro.test(a)$p.value < 0.05, "NOT NORMAL", "NORMAL")}
dact <- scan("data/dact.txt")
Normality(dact) # asmisc.r

En consecuencia, debemos aplicar para dact solo aquellos análisis y características que sean robustos a la no normalidad:

Código\(\PageIndex{11}\) (R):

dact <- scan("data/dact.txt")
summary(dact)[-4] # no mean
IQR(dact)
mad(dact)

Intervalo de confianza para la mediana:

Código\(\PageIndex{12}\) (R):

dact <- scan("data/dact.txt")
wilcox.test(dact, conf.int=TRUE)$conf.int

(Usando la idea de que cada salida de prueba es una lista, extraemos el intervalo de confianza de la salida directamente. Por supuesto, sabíamos de antemano que el nombre de un componente que necesitamos es conf.int; este conocimiento podría obtenerse de la función help (sección “Valor”). El intervalo resultante es amplio).

Para trazar datos numéricos simples, es preferible el histograma (Figura\(\PageIndex{3}\)) (las gráficas de caja son mejores para la comparación entre variables):

Figura\(\PageIndex{3}\) Histograma con curva de distribución normal superpuesta para datos dact.

Código\(\PageIndex{13}\) (R):

dact <- scan("data/dact.txt")
Histr(dact, xlab="", main="") # asmisc.r

Similar al histograma es la gráfica de vapor y hoja:

Código\(\PageIndex{14}\) (R):

dact <- scan("data/dact.txt")
stem(dact)

Además, aquí calcularemos asimetría y curtosis, tercer y cuarto momentos centrales (Figura\(\PageIndex{4}\)). La asimetría es una medida de cuán asimétrica es la distribución, la curtosis es una medida de lo puntiforme que es. La distribución normal tiene asimetría y curtosis cero mientras que la distribución uniforme “plana” tiene asimetría cero y curtosis aproximadamente\(-1.2\) (compruébalo tú mismo).

Figura Momentos\(\PageIndex{4}\) centrales (de izquierda a derecha, de arriba a abajo): por defecto, diferente escala, diferente asimetría, diferente curtosis.

¿Qué pasa con los datos dact? A partir del histograma (Figura\(\PageIndex{3}\)) y tallo y hoja podemos predecir asimetría positiva (asimetría de distribución) y curtosis negativa (distribución más plana de lo normal). Para verificar, primero es necesario cargar la biblioteca e1071:

Código\(\PageIndex{15}\) (R):

dact <- scan("data/dact.txt")
library(e1071)
skewness(dact)
kurtosis(dact)

Respuesta a la pregunta sobre los nenúfares. Primero, necesitamos verificar los datos, cargarlos en R y verificar el objeto resultante:

Código\(\PageIndex{16}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
str(ny) # asmisc.r

(La función Str () muestra los números de columna y la presencia de NA.)

Una de las posibles formas de proceder es examinar las diferencias entre especies por cada carácter, con cuatro parcelas de caja pareadas. Para hacerlos en una fila, vamos a emplear para () ciclo:

Código\(\PageIndex{17}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
oldpar <- par(mfrow=c(2, 2))
for (i in 2:5) boxplot(ny[, i] ~ ny[, 1], main=names(ny)[i])
par(oldpar)

(Aquí no, sino en muchos otros casos, para () en R es mejor reemplazar con comandos de la familia apply (). La función Boxplot acepta argumentos “ordinarios” pero en este caso, la interfaz de fórmula con tilde es mucho más práctica.)

Por favor revise esta parcela usted mismo.

Sin embargo, es aún mejor comparar caracteres escalados en una trama. La primera variante es cargar la biblioteca de celosía y crear una gráfica de enrejado similar a la Figura 7.1.8 o la Figura 7.1.7:

Código\(\PageIndex{18}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
library(lattice)
ny.s <- stack(as.data.frame(scale(ny[ ,2:5])))
ny.s$SPECIES <- ny$SPECIES
bwplot(SPECIES ~ values | ind, ny.s, xlab="")

(Como de costumbre, las parcelas de enrejado “quieren” una interfaz de forma larga y fórmula).

Por favor, revisa esta parcela tú mismo.

Alternativa es el comando Boxplots () (Figura\(\PageIndex{5}\)). No es una trama enrejada, sino diseñada con un objetivo similar de comparar muchas cosas a la vez:

Código\(\PageIndex{19}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
boxplot(ny[, 2:5], ny[, 1], srt=0, adj=c(.5, 1)) # asmisc.r

(Por defecto, Boxplots () rota las etiquetas de caracteres, pero este comportamiento no es necesario con 4 caracteres. Esta gráfica usa scale () por lo que el eje y, por defecto, no se proporciona.)

Figura Gráficas de caja\(\PageIndex{5}\) agrupadas con la función Boxplots ().

O, con un Linechart aún más nítido () (Figura\(\PageIndex{6}\)):

Código\(\PageIndex{20}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
Linechart(ny[, 2:5], ny[, 1], se.lwd=2) # asmisc.r

(A veces, los IQR son mejores para percept si agregas grid () a la trama. Pruébalo tú mismo.)

Figura Medianas\(\PageIndex{6}\) agrupadas e IQR con función Linechart ().

Evidentemente (después de SEPALOS), PÉTALOS y ESTAMENTOS hacen la mejor resolución de especies. Para obtener valores numéricos, es mejor verificar primero la normalidad.

Tenga en cuenta que la identidad de especies es la característica natural e interna de nuestros datos. Por lo tanto, es teóricamente posible que el mismo carácter en una especie exhiba una distribución normal mientras que en otra no. Es por ello que se debe verificar la normalidad por carácter por especie. Esta idea se acerca al concepto de efectos fijos que son tan útiles en los modelos lineales (ver capítulos siguientes). Los efectos fijos se oponen a los efectos aleatorios que no son naturales para los objetos estudiados (por ejemplo, si se muestrea solo una especie de nenúfares en el lago dos veces).

Código\(\PageIndex{21}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
Normality <- function(a) {ifelse(shapiro.test(a)$p.value < 0.05, "NOT NORMAL", "NORMAL")}
aggregate(ny[, 3:4], by=list(SPECIES=ny[, 1]), Normality) # asmisc.r

(Function aggregate () no sólo aplica función anónima a todos los elementos de su argumento, sino que también la divide sobre la marcha con por lista de factor (es). Similar es tapply () pero funciona sólo con un vector. Otra variante es usar split () y luego apply () la función de informes a cada parte por separado.)

Por cierto, el código anterior es bueno para aprender pero en nuestro caso particular, ¡no se requiere verificación de normalidad! Esto se debe a que los números de pétalos y estambres son caracteres discretos y por lo tanto deben tratarse con métodos no paramétricos por definición.

Así, para los intervalos de confianza, debemos proceder con métodos no paramétricos:

Código\(\PageIndex{22}\) (R):

ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")
aggregate(ny[, 3:4],by=list(SPECIES=ny[, 1]),function(.x) wilcox.test(.x, conf.int=TRUE)$conf.int)

Los intervalos de confianza reflejan la posible ubicación del valor central (aquí la mediana). Pero aún necesitamos reportar nuestros centros y rangos (¡el intervalo de confianza no es un rango!). Podemos usar tanto summary () (pruébalo tú mismo), o algún resultado personalizado que, por ejemplo, puede emplear la desviación absoluta mediana:

Código\(\PageIndex{23}\) (R):

aggregate(ny[, 3:4], by=list(SPECIES=ny[, 1]), function(.x) paste(median(.x, na.rm=TRUE), mad(.x, na.rm=TRUE), sep="±"))
ny <- read.table("http://ashipunov.info/shipunov/open/nymphaeaceae.txt", h=TRUE, sep="	")

Ahora podemos dar la respuesta como “si hay 12—16 pétalos y 100—120 estambres, esto es probablemente un lirio de agua amarillo, de lo contrario, si hay 23-29 pétalos y 66—88 estambres, esto es probablemente un lirio de agua blanco”.

Respuesta a la pregunta sobre filotaxis. Primero, necesitamos mirar el archivo de datos, ya sea con url.show (), o en la ventana del navegador y determinar su estructura. Hay cuatro columnas separadas por tabuladores con encabezados, y al menos la segunda columna contiene espacios. En consecuencia, necesitamos decirle a read.table () tanto el separador como los encabezados y luego verificar inmediatamente la “anatomía” del nuevo objeto:

Código\(\PageIndex{24}\) (R):

phx <- read.table("http://ashipunov.info/shipunov/open/phyllotaxis.txt", h=TRUE, sep="	")
str(phx)

Como ve, tenemos 11 familias y por lo tanto 11 proporciones para crear y analizar:

Código\(\PageIndex{25}\) (R):

phx <- read.table("http://ashipunov.info/shipunov/open/phyllotaxis.txt", h=TRUE, sep="	")
phx10 <- sapply(1:10, Phyllotaxis)
phx.all <- paste(phx$N.CIRCLES, phx$N.LEAVES, sep="/")
phx.tbl <- table(phx$FAMILY, phx.all %in% phx10)
dotchart(sort(phx.tbl[,"FALSE"]/(rowSums(phx.tbl))),lcolor=1, pch=19)

Figura\(\PageIndex{7}\) Dotchart muestra proporciones de fórmulas no clásicas de filotaxis.

Aquí creamos 10 primeras fórmulas clásicas de filotaxis (diez es suficiente ya que las fórmulas de orden superior son extremadamente raras), luego hicimos estas fórmulas (clásicas y no clásicas) a partir de datos y finalmente hicimos una tabla a partir de la expresión lógica que verifica si las fórmulas del mundo real están presentes en el clásico hecho artificialmente secuencia. Dotchart (Figura\(\PageIndex{7}\)) es probablemente la mejor manera de visualizar esta tabla. Evidentemente, Onagraceae (familia onagra) tiene la mayor proporción de FALSE, ahora necesitamos proporciones reales y finalmente, prueba de proporción:

Código\(\PageIndex{26}\) (R):

phx <- read.table("http://ashipunov.info/shipunov/open/phyllotaxis.txt", h=TRUE, sep="	")
phx.all <- paste(phx$N.CIRCLES, phx$N.LEAVES, sep="/")
phx10 <- sapply(1:10, Phyllotaxis)
phx.tbl <- table(phx$FAMILY, phx.all %in% phx10)
mean.phx.prop <- sum(phx.tbl[, 1])/sum(phx.tbl)
prop.test(phx.tbl["Onagraceae", 1], sum(phx.tbl["Onagraceae", ]), mean.phx.prop)

Como ves, la proporción de fórmulas no clásicas en Onagraceae (casi 77%) es estadísticamente diferente de la proporción promedio de 27%.

Respuesta a la pregunta de la encuesta de salida del “Prólogo”. Aquí está la manera de calcular cuántas personas podríamos querer pedir para estar seguros de que nuestra muestra 48% y 52% son “reales” (representan la población):

Código\(\PageIndex{27}\) (R):

power.prop.test(p1=0.48, p2=0.52, power=0.8)

¡Tenemos que preguntar a casi 5 mil personas!

Para calcular esto, se utilizó una especie de prueba de potencia que se utilizan con frecuencia para planificar experimentos. Hicimos power=0.8 ya que es el valor típico del poder utilizado en las ciencias sociales. El siguiente capítulo da definición de potencia (como término estadístico) y algo más de información sobre la salida de prueba de potencia.