Primitivas en GNU R

$ R                           # Para entrar en el entorno R una vez instalado
							  # es necesario ejecutar la R desde la consola
getwd()                       # Devuelve el
											  # directorio de trabajo
setwd("/home/davidam/public_html/software/R") # Cambia el directorio
											  # de trabajo (debe
											  # existir previamente)
list.files()                      # Lista el contenido
											  # del directorio
											  # de trabajo. Tambien ls()

source("/home/davidam/public_html/software/R/ejercicios-descriptivos.R")

q()                                           # Salir del entorno

help (aov)
?mean
help.search('mean')
??'standard deviation'
help.start()             # abre una ventana del navegador
example (mean)           # ejemplos del uso de mean

5+7      # se evalua y el resultado se imprime
			 # por pantalla pero no se guarda
5**7
5-7
5/7
5^7          # elevado a
51 %% 7  # modulo (resto de una division entera)
51 %/% 7 # division entera

x <- 1   # asigna el valor 1 a la variable x
x        # lo escribe por pantalla
y <- x^2
x = 1    # una alternativa
1 -> x   # otra alternativa

q    # Muestra informacion de la funcion
?q
# Algunas funciones elementales: log(), log10(),
# exp(), sqrt(), sin(), cos(), tan()
# Mas funciones basicas: c(), max(), min(),
# pmax(), pmin(), range(), length(), sort()
# sum(), rowSums(), prod(), mean(), colMeans(),
# var(), cumsum(), etc.
c(2,69,5,47,3647)
sum(1,3,5,7)
raiz64 <- sqrt(64)   # el resultado de una funcion
					 # puede guardarse en una variable
print(x)
print(pi,digits=16); print(pi,digits=5)

ls()
rm(x,raiz64)
ls()
rm(list=ls())
ls()

verdad <- T
mentira <- F
edad <- c(27,30,29);
is.numeric(edad)
is.integer(edad)
is.double(edad) # tambien existen is.complex(),
				# is.logical(), is.character()...
typeof(edad)    # devuelve el tipo del objeto
				# introducido como argumento

calle <- c("c/ rue del percebe","28013");

as.character(verdad)
as.numeric(verdad)      # doble precision por defecto
y <- c(0,1,27,14.9,NA,-3.2);y   # ojo, trunca en vez de redondear
as.logical(y)
as.character(y)
as.integer(y)
y<-c("hola", "adios","true","TRUE", "T", "F","129",NA ,"6");y
as.logical(y)
as.numeric(y)

x <- c(1,4,6,3,7); x
x[3]            # devuelve el valor guardado en 3a posicion
x[-4]           # devuelve todo excepto el valor
						# guardado en 4a posicion
x[1] <- 35      # guarda el valor 35 en la primer posicion
length(x)       # devuelve la longitud de x
x[11] <- 99; x

x <- c(1,2,3,4,5)
y <- c(5,4,3,2,1)
x+y
y <- c(1,2,3);y
suma <- x+y; suma       # recicla el vector mas corto
suma2 <- x-1            # con escalares parece mas natural
x <- c(0.0001, 0.001, 0.01, 0.1, 1)
log10(x)            # muchas funciones tambien se
								# aplican vectorizadamente

1:10
1:10-1              # el operador ':' tiene maxima prioridad
1:{10-1}
#funcion seq(from, to, by, length)
x <- seq (-3,3,0.5);x       # valores desde -3 hasta 3 a intervalos de 0.5
seq (-3,3,length=10)        # 10 valores equiespaciados desde -3 hasta 3
seq (-3, by=0.5,length=10)  # 10 valores a intervalos de 0.5 desde -3
seq(along=x)            # Genera la secuencia 1, 2, 3,..., length(x)
#funcion rep(x, times, each)
x <- 1:4
rep(x,times=3)          # repite el contenido de x tres veces
rep(x,each=3)               # repite cada elemento de x tres veces
rep(x,times=2,each=3)
rep(c(0,1),times=c(4,3))    # 'times' puede ser un vector. Aqui ya no cabe 'each'
x <- seq(19,6,-3);x
rep(x,1:length(x))      # 'times' puede estar implicito times={1:length(x)}

x <- c(20,80,30,50,0)
order (x, decreasing=F)     # devuelve las posiciones del vector ordenadas segun su contenido
sort (x, decreasing=F)      # devuelve el contenido del vector ordenado
rank(x)             # devuelve el orden de cada posicion segun su contenido
min(x)
which.min(x)            # equivale a which(x == min(x))
x <- c(1,1,3:1,1:4,3); y <- c(0,9:1)
x.ord <- order(x,y)     # ordena los valores de x y, en caso de empate, utiliza los valores de y en las posiciones correspondientes. Ambos, x e y deben tener la misma longitud
x.ord;x;y

x <- c(1:10); x
valor.verdad <- x>5; valor.verdad
x <- 1:10
x[x >= 5] <- 20; x  # condicional implicito
x[x == 1] <- 25; x
x <- 1:3; y <- 1:3
x == y          # equivalente a identical(x, y) o all.equal(x, y, tol=0)
y <- y+0.001; y
identical(x, y)
all.equal(x, y, tol=0.001)

a<-c(TRUE,TRUE,FALSE,FALSE)
b<-c(F, T, F, T)
a&b
a|b
!(a&b)
(!a)|(!b)

rm(list=ls())
x <- c(3,6,4,2,8,9)
print (x); length(x)
x[8:10] <- 3;x
is.na(x)
!is.na(x)
which(is.na(x))
x[is.na(x)]<-999;x  # codifica como 999 los valores perdidos
x==NA           # la expresion logica x == NA sa un resultado muy distinto de
is.na(x)

x <- c(0,7,8); x/x
1/x
-1/x
is.nan(x)
is.nan(x/x)
is.nan(1/x)
# Hay que tener cuidado con los NaN porque cualquier operacion con un NaN resulta en un NaN

# Concatena objetos en un vector de caracteres. Funcion paste(..., sep = " ", collapse = NULL)
juntar <- paste("Una ", "frase ", "cualquiera",  collapse ="");juntar
v1<-c("A","B")
v2<-2:3
codigos <- paste(v1,v2, sep = "");print(codigos)
codigos <- paste(v1,v2, sep = ".");print(codigos)
x <- paste(LETTERS[1:5]);x
x <- paste(LETTERS[1:5], collapse="");x

# Concatena e imprime. Funcion cat(... , file = "", sep = " ", fill = F, labels = NULL, append = F)
verano <- month.abb[7:9]; verano
cat(verano)         # el resultado no puede guardarse en ua variable
cat(verano,"\n")
cat(verano, sep=',')        # concatena separando con comas e imprime
cat(verano, sep=';', fill=3)
cat("Estaciones:","\t","Moncloa","\n", "\t", "\t","Aluche","\n")

#Todo junto: cat() y paste()
x<-2/3; cat(paste("resultado", signif(x,2), sep=" : " ),"\n")

edad <- c(12,22,15,16,10)
names(edad)                     # por defecto no se asignan nombres
names(edad) <- paste("suj#", sep = "",length(edad):1); edad
names(edad)
edad["suj#2"]                   # devuelve la edad de Sujeto #2 (almacenado en la penultima posicion)
edad[length(edad)-1]
names(edad) <- NULL; edad   # elimina los nombre asignados

X <- 1:12
dim(X) <- c(4,3); X # los elementos se organizan por columnas
matrix(1:12,nrow=3,ncol=4,byrow=T)
matrix(1:12,nrow=4,byrow=T)
Mi.matriz <- matrix(1:12,3,4,F); Mi.matriz
tamano <- dim(Mi.matriz)     # otro uso de dim(); devuelve las dimensiones de la matriz
rownames(Mi.matriz) <- LETTERS[1:tamano[1]]; Mi.matriz
colnames(Mi.matriz) <- paste("Var",1:tamano[2], sep=""); Mi.matriz
dimnames(Mi.matriz)

X1 <- c(3,7,5)
X2 <- c(8,3,1)
Xx <- cbind(X1,X2);Xx
Yy <- rbind(X2,X1);Yy
Zz <- cbind(X1,Xx);Zz

X <- matrix(c(1,4,12,15),2,2);X
X[1,2]               # elemento guardado en la 1a fila, 2da columna
X[1, ]               # todos los elementos de la primera fila
X[ ,2]               # todos los elementos de la segunda columna
X[3]                 # para pensar un poco...
Mi.matriz['B',]  #tambien se pueden usar los nombres (se han asignado)

X <- matrix(c(1,4,12,15),2,2); X
Y <- matrix(1:4,2,2); Y
X+Y
X-Y
X%*%Y                        # producto matricial
X*Y                      # producto elemento por elemento
t(X)                     # traspuesta
det(X)                   # determinante
X.inv <- solve (X)   # inversa de X (siempre que X sea cuadrada no singular, claro)
X.inv
X%*%X.inv                # Comprueba el resultado. Ojo a los errores de redondeo
# En general, solve(a,b) es una funcion que resuelve a %*% x = b para x, donde b puede ser un vector o una matriz. Si no se explicta se asume que es la matriz identidad y la funcion devuelve la inversa de a
A <- matrix(c(1,4,12,15),2,2); A
B <- matrix(c(5,2),2,1);B
X <- solve(A,B); X
A%*%X            # comprueba que el resultado es efectivamente B

X <- matrix(c(1,4,12,15),2,2); X
rowSums(X)
colSums(X)
rowMeans(X)
colMeans(X)

apply(X,1,sum)       # suma por filas; equivale a rowSums(X)
apply(X,2,mean)  # media por columnas; equivale a colMeans(X)
apply(X,1,sd)

A <- array(1:24, c(3, 4, 2)); A
dimnames(A) <- list(c("fila1", "fila2", "fila3"), c("col1", "col2", "col3", "col4"), c("capa 1", "capa 2"));A

length(lista1)       # devuelve el numero de componentes de la lista
str(lista1)      # devuelve informacion sobre la estructura de la lista
a <- lista1[[1]] # devuelve el 1er objeto de la lista en forma de vector/matriz (y nombre  excluido)
a; typeof(a)
b <- lista1[1]       # devuelve una sublista compuesta por los elementos de la 1a entrada de la lista (nombre incluido)
b; typeof(b)

nomb <- c("reales","enteros", "matriz", "caracteres", "logicos")
names(lista1) <- nomb;   # asigna nombres a los componenes de la lista
names(lista1)                # devuelve los nombres de los componentes (si los hay)
lista1[["reales"]]
lista1$logicos           # equivalente a lista1[["logicos"]] y a lista[[5]]

# Otro modo de definir una lista que incluye nombres para los componentes
lista2 <- list(A=mis.num, B=mis.num2); lista2
names(lista2)

attach(lista1, warn.conflicts = T)

rownames(mis.datos); colnames(mis.datos)
dim(mis.datos)               # dimensiones
str(mis.datos)               # estructura
attributes(mis.datos)    # clase y nombres de fila y columna

trozos <- split (mis.datos, mis.datos$var5); trozos      # devuelve una lista "trozos" con los valores de var1 segregados segun var5
typeof(trozos)
junto.otravez <- unsplit(trozos, mis.datos$var5); junto.otravez

# el 2do argumento funciona como factor, no es necesario que este definido como tal.
# Ademas se pueden seleccionar solo algunas variables (ej, var2 y var3)
trozos <- split (c(mis.datos$var1,mis.datos$var2), mis.datos$var3); trozos

# lo mismo podria conseguirse con una sucesion de instrucciones como esta:
twitter.si <- mis.datos[mis.datos$var5=="UsaTwitter",];twitter.si

subset(mis.datos, var3 >= 15)
subset(mis.datos, var3 >= 15, select=c(var1,var2))
subset(mis.datos, var3 >= 15, select=-var3)

subset(data$q_005, (data$q_005 > 1900) & (data$q_005 < 2013))

mis.datos2 <- transform (mis.datos, log.var1= log(var1)); mis.datos2

edit(data.frame())            # abre el editor con una 'hoja de datos' vacia
					  # los huecos se rellenaran solos con NA
					  # doble click en "var?" para cambiar el nombre de la variable
					  # ¿que ha pasado?, ¿donde estan los datos?
 respuestas <- edit(data.frame()) # Es esencial asignarle un nombre a la hoja de datos
 respuestas
 respuestas2 <- edit(respuestas)  # guarda los cambios en una nueva hoja de datos (y conserva la antigua)
 respuestas;respuestas2
 fix(respuestas)                      # sobreescribe los cambios. Es equivalente a respuestas <- edit(respuestas)

search()     # muestra la ruta de busqueda
# Cuidado con los nombres de las variables del espacio de trabajo: Tienen prioridad
rm(list=ls())
women
attach(women)
mean(weight)
weight <- 0; mean(weight)
mean(women$weight)

rm(list=ls())
data(women)
write.table(women, "../datos/mujeres.dat", append=T) # Si la carpeta datos no existe previamente dara un mensaje de error
?write.table
write.table(women, "../datos/mujeres.dat", append=T, row.names=F, col.names= F, sep='', dec=',')
write.table(women, "../datos/mujeres.dat", append=F,quote=F)

x<-seq(1,3,0.1)
write(x,"../datos/numeros.dat", ncolumns=1,append=T)
write(x,"../datos/numeros.dat", ncolumns=length(x),append=F)

?read.table
rm(list=ls())
mis.datos <- read.table("../datos/mujeres.dat", header=T, check.names=T, sep=" ",na.strings=999, comment.char = "#")
edit(mis.datos)
median(mis.datos$weight)
colMeans(mis.datos)
attach(mis.datos)
sum(weight >=125)/length(weight)     # calcula la proporcion de individuos que pesan al menos 125 lb

rm(list=ls(all=TRUE))
b=read.table("Oster.txt", header=T)
dim(b)
summary(b)
names(b)# obtengo los nombres de las cabeceras de columna
head(b)# Quiero ver los 6 primeros valores
tail(b) # Los seis últimos
attach(b)# anclamos la base de datos para trabajar con ella
View(b)# Veo mis datos
edit(b)

library(foreign) # En primer lugar hay que cargar el paquete
mi.encuesta.spss <- read.spss("../datos/SurveyStatisticsI.sav", to.data.frame=T)
edit(mi.encuesta.spss)

?scan
mis.datos <- scan("../datos/mujeres.dat", what = list(character(), numeric(),numeric()), skip=1);mis.datos
mis.datos2 <- scan("../datos/mujeres.dat", what = list(0,0,0), skip=3)

x11()                                    # abre una nueva ventana grafica (Unix y Wndows)
windows(width=15, height=15, title='Figura 1')   # abre una nueva ventana grafica (Windows)
quartz()                             # abre una nueva ventana grafica (Mac OS)
pdf("prueba.pdf")                    # abre un pdf de nombre "prueba.pdf" y redirige el output gráfico al archivo pdf
postscript("prueba.eps")     # abre un archivo de postscript encapsulado de nombre "prueba.eps" y redirige el output gráfico

dev.list()   # lista todos los dispositivos grafico abiertos. Los nums son los identificadores
dev.cur()    # devuelve el identificador del dispositivo grafico activo
dev.set(2)   # establece como dispositivo gráfico activo el dispositivo num 2
dev.off()    # cierra el dispositivo grafico activo; dev.off(3) cierra el dispositivo grafico num 3

x <- 1:15; y <- 30:16
# con datos cuantitativos produce diagramas de dispersion
plot (x,y, pch=2)
plot(women)      # puede usarse una hoja de datos o matriz. OJO solo toma las 2 primeras columnas
plot(y)          # representa los valores de 'y' en funcion de su posicion en el vector

# plot (factor) produce un diagrama de barras
f <- factor(rep(c(1,2,2,1,1,3), each=2), levels=1:3, labels=c("No fumador", "Fumador", "NS/NC"))
plot(f)
# tambien existe la funcion barplot para crear diagramas de barras

# plot (var.cuant, factor) produce diagramas de caja organizados según los niveles del factor
entorno <- factor(rep(c(1,2,2), each=5), levels=1:2, labels=c("Urbano", "Rural"))
plot(entorno,y)

x <- 1:15; y <- 30:16
# con datos cuantitativos produce diagramas de dispersion
png(filename="/tmp/figure.png", height=295, width=300, bg="white")
plot (x,y, pch=2)
dev.off()

datos2 <- data.frame(entorno,women)
coplot(women$weight ~ women$height | entorno) # representa peso en funcion de altura separando segun entorno
socioeconomico<- factor(c(3,3,2,2,2,2,1,1,1,1,1,1,1,2,2), levels=1:3, labels=c("bajo","medio","alto"))
plot(socioeconomico)
coplot(women$weight ~ women$height | entorno+socioeconomico)

pers <- matrix(rnorm(1000), ncol = 5)    # genera 1000 datos con distribucion normal. Vease seccion 'Muestreo y probabilidad'
colnames(pers) <- c("extraversion", "estabilidad", "apertura", "responsabilidad", "amabilidad")
pairs(pers)
# Nótese la diferencia con el resultado de plot()
x11(); plot(pers)            # usa solo las 2 primeras variables

pers <- as.data.frame(matrix(rnorm(1000), ncol = 5))
colnames(pers) <- c("extraversion", "estabilidad", "apertura", "responsabilidad", "amabilidad")
cor(pers, method='pearson')  # matriz de correlaciones. Otros metodos: 'spearman' y 'kendall'
round(cor(pers),3)   # la misma matriz de antes, pero con solo 3 decmales
attach(pers)
hist(extraversion)
x11(); hist(extraversion, labels=T, col= "gray", xlim=c(-4,4), ylab="Frecuencia absoluta", sub= "Figura 1. Ejemplo de histograma", main= "Extraversión")
hist(amabilidad, col= "red", border="white", add=T, xlim=c(-4,4), sub= "Figura 2. Otro ejemplo de histograma", main= c("Extraversión  Amabilidad superpuestos"))
x11(); hist(extraversion,freq=F, labels=T, col= "light blue",xlim=c(-4,4),ylab="Frecuencia relativa", las=1,sub= "Figura 1b. Ejemplo de histograma", main= "Extraversión")
hist(amabilidad, col= "gray", border="white", xlim=c(-4,4), sub="", main= "Figura 3. Otro ejemplo de histograma cambiando el color")

pdf("ejemplos_histograma.pdf")
hist(extraversion, nclas=3, col= "light blue", sub= "Figura 3. Ejemplo de histograma", main= "Extraversión")
hist(extraversion, breaks=seq(-4,4,0.5),col= "pink",main= "Figura 3. Ejemplo de histograma")
dev.off()

boxplot(extraversion, main= "Diagrama de caja de extraversión")
text(locator(5),"más outliers", adj=0)
boxplot(pers)

qqnorm(extraversion, col="red")
qqline(extraversion)

colors()

   ?par
   par.actual <- par()      # guarda los valores actuales (para poder reestablecerlos posteriormente)
   par(font.main=4, las=1, xaxs= "r", col.lab="blue", col="skyblue3",lwd=3)
   hist(extraversion,xlim=c(-4,4))
   plot(extraversion, amabilidad)
#   par(par.actual)         # reestablece los valores guardados
   par(mfrow=c(2,1))        # distribuye el espacio grafico e dos filas

curve(x^2)
curve(x^2, from=-3, to=3) # los cambios introducidos a traves de par no tienen efecto en la nueva ventana
# tambien podria haberse hecho con plot
x11()
x <- seq(-3,3,0.25)
y <- x^2
plot(x,y)
plot (x,y, type="o") # 'p': puntos, "l": lines, 'b': puntos conectados mediante lineas, etc.
# un ultimo ejemplo
hist(extraversion, freq=F)
curve(dnorm(x), add=T)
detach(pers)

datos <- c(1,5,9,8,52,2,4)
sample(datos,3)              # toma 3 elementos de 'datos'. Por defecto, muestrea SIN reposición

dado <- 1:6                      # para simular el lanzamiento de 5 dados hace falta muestreo CON reposición
sample(dado,5,replace=T)
# sample() tambien permite incluir un vector de pesos si los sucesos no son equiprobables

baraja <- paste(rep(c("As",2:7,"Sota","Caballo","Rey"),4),rep(c("Oros","Copas","Espadas","Bastos"),each=10));baraja
sample(baraja,2)

#  set.seed() # inicializa el generador de aleatorios con una semilla concreta (reproductibilidad)
   set.seed(-98)

datos.gauss <- rnorm (100, mean=5, sd=2) # genera 100 numeros ~N(5,2)
runif (4,-10,10)                                 # genera 4 numeros ~U(-10, 10)

x1 <- seq(-5,5,0.01)
gauss <-dnorm(x1)
plot (x1,gauss,"l", lwd=2, col="blue")
student <- dt(x1, df=4)
x11();
plot (x1,student,"l", lwd=2,col="green")
x2 <- seq(0,10,0.01)
chi2 <-dchisq(x2,df=3)
x11();
plot (x2,chi2,"l",lwd=2,col="red")

pnorm(0.7)                               # devuelve la prob de que una V.A.con distribución N(0,1) tome un valor menor o igual que 0.7
pnorm(0.7, mean=10, sd=0.5)
pchisq(3.84, df=3, lower.tail=F) # devuelve la probabilidad de que una V.A.~ Chi2 con 3gl tome un valor mayor que 3.84

qnorm(0.025)                      # devuelve el cuantil 2.5 (valor de la V.A. que deja por debajo el 2.5% de la distribución)
qf(0.975, df1=4, df2=53)      # se puede usar para obtener los puntos de corte de la región de rechazo en un contraste de hipótesis

> dbinom(4, size=12, prob=0.2)
[1] 0.1329
> pbinom(4, size=12, prob=0.2)
[1] 0.92744

> ppois(16, lambda=12)   # lower tail
[1] 0.89871
> ppois(16, lambda=12, lower=FALSE)   # upper tail
[1] 0.10129

> runif(10, min=1, max=3)
[1] 1.6121 1.2028 1.9306 2.4233 1.6874 1.1502 2.7068
[8] 1.4455 2.4122 2.2171

> pexp(2, rate=1/3)
[1] 0.48658

> pnorm(84, mean=72, sd=15.2, lower.tail=FALSE)
[1] 0.21492

> qchisq(.95, df=7)        # 7 degrees of freedom
[1] 14.067

> qt(c(.025, .975), df=5)   # 5 degrees of freedom
[1] -2.5706  2.5706

> qf(.95, df1=5, df2=2)
[1] 19.296

> ppois(16, lambda=12)   # lower tail
[1] 0.89871

> ppois(16, lambda=12, lower=FALSE)   # upper tail
[1] 0.10129

binom.test(51, 336, (1/6))

scipy.stats.binom_test(51, 336, 1.0/6)

> qf(.95, df1=5, df2=2)
[1] 19.296

> library(MASS)
> levels(survey$Smoke)
> smoke.freq = table(survey$Smoke)
> smoke.freq
> smoke.prob = c(.045, .795, .085, .075)
> chisq.test(smoke.freq, p=smoke.prob)

		Chi-squared test for given probabilities

data:  smoke.freq
X-squared = 0.1074, df = 3, p-value = 0.991

> library(MASS)       # load the MASS package
> tbl = table(survey$Smoke, survey$Exer)
> tbl                 # the contingency table

		Freq None Some
  Heavy    7    1    3
  Never   87   18   84
  Occas   12    3    4
  Regul    9    1    7

> chisq.test(tbl)

		Pearson’s Chi-squared test

data:  table(survey$Smoke, survey$Exer)
X-squared = 5.4885, df = 6, p-value = 0.4828

Warning message:
In chisq.test(table(survey$Smoke, survey$Exer)) :
  Chi-squared approximation may be incorrect

> library(MASS)         # load the MASS package
> head(immer)
  Loc Var    Y1    Y2
1  UF   M  81.0  80.7
2  UF   S 105.4  82.3

> wilcox.test(immer$Y1, immer$Y2, paired=TRUE)

		Wilcoxon signed rank test with continuity correction

data:  immer$Y1 and immer$Y2
V = 368.5, p-value = 0.005318
alternative hypothesis: true location shift is not equal to 0

Warning message:
In wilcox.test.default(immer$Y1, immer$Y2, paired = TRUE) :
  cannot compute exact p-value with ties

> am.glm = glm(formula=am ~ hp + wt, data=mtcars, family=binomial)
> newdata = data.frame(hp=120, wt=2.8)
> predict(am.glm, newdata, type="response")

> am.glm = glm(formula=am ~ hp + wt,
+              data=mtcars,
+              family=binomial)

> summary(am.glm)

Call:
glm(formula = am ~ hp + wt, family = binomial, data = mtcars)

Deviance Residuals:
	Min       1Q   Median       3Q      Max
-2.2537  -0.1568  -0.0168   0.1543   1.3449

Coefficients:
			Estimate Std. Error z value Pr(>|z|)
(Intercept)  18.8663     7.4436    2.53   0.0113 *
hp            0.0363     0.0177    2.04   0.0409 *
wt           -8.0835     3.0687   -2.63   0.0084 **
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

(Dispersion parameter for binomial family taken to be 1)

	Null deviance: 43.230  on 31  degrees of freedom
Residual deviance: 10.059  on 29  degrees of freedom
AIC: 16.06

Number of Fisher Scoring iterations: 8

boxplot(Peso ~ Manejo)
wilcox.test(Peso~Manejo)
boxplot(Ostertagia ~ Manejo)
plotmeans(Ostertagia~Manejo,col="coral", connect=T, ylab="Carga parasitaria (Ostertagia)", xlab="Manejo", data=b)# plot means and CI95% Grafica mitjana

> df1 = read.table("fastfood-1.txt", header=TRUE); df1
  Item1 Item2 Item3
1    22    52    16
2    42    33    24
3    44     8    19
4    52    47    18
5    45    43    34
6    37    32    36

> r = c(t(as.matrix(df1))) # response data
> r
 [1] 22 52 16 42 33 ...

> f = c("Item1", "Item2", "Item3")   # treatment levels
> k = 3                    # number of treatment levels
> n = 6                    # observations per treatment

> tm = gl(k, 1, n*k, factor(f))   # matching treatments
> tm
 [1] Item1 Item2 Item3 Item1 Item2 ...

> av = aov(r ~ tm)
> summary(av)
			Df Sum Sq Mean Sq F value Pr(>F)
tm           2    745     373    2.54   0.11
Residuals   15   2200     147