# VERSI_N EN TEXTO DEL C_DIGO EN R PARA LA CLASE
# DE SIGNIFICANCIA EN R. OCTAVIO MART_NEZ MARZO 2011

# Solamente se incluye lo que se debe teclear,
# no los resultados que da R.

A <- c(587941, 89933, 9160694-587941, 5743943-89933)
# vector que contiene los cuatro n£meros que necesitamos

A <- matrix(A, nrow=2, ncol=2, byrow=T, dimnames = list(c("A","NoA"), c("To","In")))


# Experimentando com matrices
matrix(c(1,2,3,4), nrow=2, ncol=2)
matrix(c(1,2,3,4), nrow=2, ncol=2, byrow="T")

# Experimentando con listas
list( c("A","NoA"), c("To","In"))

# Sumas de A
sum(A)
apply(A,1,sum) # "sum" se aplica a la primera dimensi¢n (renglones)
apply(A,2,sum) # "sum" se aplica a la segunda dimensi¢n (columnas)
f.A <- A[1,]/apply(A,2,sum) # Frecuencia relativa de A
f.A

# Resultado "erroneo":
A/apply(A,2,sum)

rep(apply(A,2,sum), each=2) # Cada elemento del vector repetido dos veces
f.A <- A/rep(apply(A,2,sum), each=2) # Ahora si, nos da f[A|To], ... etc 
f.A
apply(f.A,2,sum)

# An lisis de los datos
f.A[1,1]
f.A[1,2]
f.A[1,1]/f.A[1,2]

pe.A <- A/sum(A) # Dividimos cada celda por el total...
pe.A
apply(pe.A,1,sum) # Probabilidades estimadas de A y NoA...
apply(pe.A,2,sum) # Probabilidades estimadas de To e In...
apply(pe.A,1,sum) %*% t(apply(pe.A,2,sum)) 
pe.A.H0 <- apply(pe.A,1,sum) %*% t(apply(pe.A,2,sum)) 
# Ahora le ponemos el atributo de nombres de renglones:
attributes(pe.A.H0)$dimnames[[1]] <- c("A","NoA")
pe.A.H0
pe.A
f.A
apply(f.A,2,sum) # Suma por columnas de f.A
rep(apply(pe.A.H0,2,sum),each=2) # Para hacer la divisi—n...
pe.A.H0 / rep(apply(pe.A.H0,2,sum),each=2)
EA.H0 <- pe.A.H0*sum(A)
EA.H0
A # Frecuencias observadas (matriz original de datos)
A-EA.H0 # La diferencia entre observados y esperados

# Densidad Chi-cuadrad
plot(c(1:1000)/100, dchisq(c(1:1000)/100,df=1), type="l", xlab="x", ylab="f(x) - Densidad Chi Cuadrada",  main="Grafica de densidad Chi Cuadrada con un grado de libertad", col="red")

qchisq(.95,df=1) # Valor para 1-alfa=0.95; alfa=0.05
qchisq(.99,df=1) # Valor para 1-alfa=0.99; alfa=0.01
qchisq(.999,df=1) # Valor para 1-alfa=0.999; alfa=0.001

qchisq(1-0.000001,df=1) # Valor "cr¡tico" de chi cuadrada

# Obtenci¢n del estad¡stico
A-EA.H0
(A-EA.H0)^2 # Nos da los valores de (Oi - Ei) al cuadrado
((A-EA.H0)^2)/EA.H0 # Nos da los valores de (Oi - Ei)^2 / Ei
sum(((A-EA.H0)^2)/EA.H0) # Nos da el valor del estad¡stico chi cuadrada

1-pchisq(191476.2, df=1) # 1 - P[Chi<=191476.2]
pchisq(191476.2, df=1, lower.tail=FALSE) # P[Chi>191476.2]

chisq.test(A) # Realizando la prueba sobre la matriz A

chi.A <- chisq.test(A) # Guardando los resultados para posterior an lisis

attributes(chi.A)$names
chi.A$expected

EA.H0

chisq.test(A, correct=FALSE) # Realizando la prueba sin la correccion de Yates.

fish.A <- fisher.test(A) # Prueba exacta de Fisher para nuestros datos

fish.A

A.dis <- A/100000 # Nuestra nueva tabla "disminuida"
A.dis # A£n cuando no tiene valores enteros (como deber¡a)

(5.87941/(5.87941+85.72753))/(0.89933/(0.89933+56.54010))

chisq.test(A.dis)
fisher.test(A.dis)


# Datos para Figura 1
factor <- 10^c(0:5) # Factor de reducci¢n de la muestra 1, 10, 100, 1000, 10000, 100000.
n <- sum(A)/factor # Aqui guardaremos el tama¤o de muestra
val.chi <- c() # Aqui guardaremos el valor del estadistico chi cuadrada
val.P <- c() # Aqui guardaremos el valor de P (significancia)

for(i in c(1:6)){ # Iniciamos un "bucle" para hacer las cuentas en cada caso.
temp <- chisq.test(A/factor[i]) # Guardamos en un objeto "temp" el resultado de la prueba
val.chi <- c(val.chi, temp$statistic) # Recuperamos el valor del estaditico
val.P <- c(val.P, temp$p.value ) # Recuperamos el valor de P 
} # Cerramos el bucle, terminando la ejecucion
varNinA <- data.frame(n,factor,val.chi,val.P) # Reunimos las columnas en un marco de datos

varNinA

# Gr fica (Figura 1)
plot(log10(varNinA$n), log10(varNinA$val.chi), ylim=c(-0.5, 5.5), xlim=c(2,7.5), xlab="Log 10 del tama¤o de muestra", ylab="Log 10 del valor del estad¡stico Chi cuadrada",main="Relaci¢n entre el tama¤o de la muestra y el valor del\nestad¡stico chi cuadrada",sub="Datos para el gen A")

points(log10(varNinA$n), log10(varNinA$val.chi), type="l", col="red")

text(log10(varNinA$n), log10(varNinA$val.chi), labels=round(varNinA$val.P,7), pos=3, cex=.7, col="blue")

text(log10(varNinA$n), log10(varNinA$val.chi), labels=round(varNinA$n), pos=1, cex=.7, col="orange")

legend("topleft", legend=c("P - significancia","n - Tama¤o de muestra"), text.col=c("blue","orange"),cex=0.7, bty="n")

# Explorando tama¤os de muestra
n <- c(100:800) # Los tama¤os de muestra que exploraremos
val.chi <- c()
val.P <- c()

for(i in c(1:701)){ # Iniciamos un "bucle" para hacer las cuentas en cada caso.
temp <- chisq.test(n[i]*pe.A) # Guardamos en un objeto "temp" el resultado de la prueba
val.chi <- c(val.chi, temp$statistic) # Recuperamos el valor del estad¡stico
val.P <- c(val.P, temp$p.value ) # Recuperamos el valor de P 
} # Cerramos el bucle, terminando la ejecuci¢n

varNinA2 <- data.frame(n,val.chi,val.P)

plot(varNinA2$n, varNinA2$val.P, type="l", col="red", xlab="n - Tama¤o de muestra", ylab="P - significancia de la prueba", main="Valores de significancia para la prueba chi cuadrada\nen la tabla de contingencia",sub="(Datos del gen A)") # Todos los valores del dataset

plot(varNinA2$n[val.P<0.01], varNinA2$val.P[val.P<0.01], type="l", col="red", xlab="n - Tama¤o de muestra", ylab="P - significancia de la prueba", main="Valores de significancia para la prueba chi cuadrada\nen la tabla de contingencia",sub="(Datos del gen A). Solamente valores de P<0.01") # Solo los valores con P<0.01


n <- 0 # Para n = 0
((587941-n)/9160694)/( (89933+n)/5743943)
[1] 4.099175
n <- 100 # Para n = 100
((587941-n)/9160694)/( (89933+n)/5743943)
[1] 4.093925
n <- 1000 # Para n = 1000
((587941-n)/9160694)/( (89933+n)/5743943)
[1] 4.0472
n <- 100000 # Para n = 100000
((587941-n)/9160694)/( (89933+n)/5743943)
[1] 1.610826

# Funci¢n para crear matriz modificada
modA <- function(n) {A + c(-n, n, n, -n)}

A # La matriz original
modA(10) # La matriz modificada con n=10

# Funci¢n para la razon o radio de probabilidades
radio <- function(X) { (X[1,1]/(X[1,1]+X[2,1])) / (X[1,2]/(X[1,2]+X[2,2])) }

radio(A)
radio(modA(1000))
radio(modA(10000))
radio(modA(100000))

EA.H0 # Valores esperados bajo independencia en nuestra matriz original

radio(EA.H0) # radio de las probabilidades esperadas 

587941 - 416635 # Cambio de A[1,1] en la observada respecto a la esperada

# Explorando cambios de razones de probabilidades
n <- 1000*c(0:172) # Valores en que cambiaremos n (de 0 a 172000)
p.ATo <- c() # Vector para guardar P[A|To]
p.AIn <- c() # Vector para guardar P[A|In]
p.dif <- c() # Vector para guardar P[A|To]-P[A|In]
p.radio <- c() # Vector para guardar P[A|To] / P[A|In]
chi.cuad <- c() # Vector para guardar el valor del estad¡stico chi cudrada
P <- c() # Vector para guardar la significancia del estad¡stico chi cudrada

for(i in 1:173) { # Para todos los valores desde i a 173 (longitud de n)
temp <- modA(n[i]) # Matriz A, pero haciendo la operaci¢n con n[i] 
p.ATo <- c(p.ATo, temp[1,1]/(temp[1,1]+temp[2,1])) # Prob estimada de A|To
p.AIn <- c(p.AIn, temp[1,2]/(temp[1,2]+temp[2,2])) # Prob estimada de A|In
p.dif <- c(p.dif, p.ATo[i]-p.AIn[i]) # Diferencias de las probabilidades
p.radio <- c(p.radio, radio(temp)) # Radio o raz¢n de las probabilidades
temp.test <- chisq.test(temp) # Objeto con los resultados de la prueba estadistica
chi.cuad <- c(chi.cuad, temp.test$statistic) # Capturamos el valor del estad¡stico
P <- c(P, temp.test$p.value) # Capturamos el valor de P
} # cerramos el bucle para ejecutar...

sign.r <- data.frame(n, p.ATo, p.AIn, p.dif, p.radio, chi.cuad, P) # Marco de datos

attach(sign.r) # Ponemos el marco de datos "a la vista"

# Graficando
plot(sign.r$p.radio, sign.r$P, type="l", xlab= "radio = P[A | To] / P[A | In]", ylab="Significancia (P)", main = "Significancia como funci¢n del radio\nen tablas con 14,904,637 etiquetas")

plot(sign.r$p.dif, sign.r$P, type="l", xlab="diferencia = P[A | To] - P[A | In]", ylab="Significancia (P)", main="Significancia como funci¢n de la diferencia de probabilidades\nen tablas con 14,904,637 etiquetas")

plot(sign.r$p.radio[sign.r$P>0], sign.r$P[sign.r$P>0], type="l",xlab="radio = P[A | To] / P[A | In]",  ylab = "Significancia (P)", main="Significancia como funci¢n del radio\nen tablas con 14,904,637 etiquetas", sub="Solo valores P>0")

plot(sign.r$p.dif[sign.r$P>0], sign.r$P[sign.r$P>0], type="l",xlab="diferencia = P[A | To] - P[A | In]", ylab="Significancia (P)", main="Significancia como funci¢n de la diferencia de probabilidades\nen tablas con 14,904,637 etiquetas", sub="Solo valores P>0")

# Constante de raz¢n de probabilidades
cons <- 5743943 / 9160694 # Que es igual a 0.6270205

tot.col.A <- apply(A,2,sum) # Los totales de las columnas de la matriz original

c(1, cons)/tot.col.A # Para b=1, a=to.in

c(10000, 10000*cons)/tot.col.A # Para a=10000,

forma <- function(a, b) {matrix(c(a, 9160694-a, b, 5743943-b), nrow=2, ncol=2,   dimnames = list(c("A", "NoA"), c("To", "In")))}

A
forma(a=587941, b=89933) # Formando una matriz id‚ntica a A
B <- round(forma(a=100000, b=100000*cons)) # Una matriz con r=1
B

chisq.test(B)

# Variando numero de etiquetas y raz¢n
prob.A.To <- 10^-c(6:2) # Las probabilidades P[A|To]

tags.a <- round(9060694 * prob.A.To) # Numeros de tags "a" correspondientes a las P[A|To]

tags.a

radios <- c(100:300)/100 # radios P[A|To]/P[A|In].
i1 <- c() # Indicador 1 (n£mero de tags)
i2 <- c() # Indicador 2 (radios)
t.a <- c() # Para guardar el n£mero de tags a
t.b <- c() # Para guardar el n£mero de tags b
tags <- c() # para la suma de etiquetas t.a + t.b
r <- c() # Para guardar el radio
p.A.To <- c() # Para guardar P[A|To]
p.A.In <- c() # Para guardar P[A|In]
chi <- c() # Para guardar el valor del estad¡stico chi
P <- c() # Para guardar P (significancia)

k <- 0
for(i in c(1:5)){ # para cada valor en tags.a
for(j in c(1:201)) { # para cada valor de "radios"
i1 <- c(i1, i) # Indicador 1 (n£mero de tags)
i2 <- c(i2, j) # Indicador 2 (radios)
k <- k+1
t.a <- c(t.a, tags.a[i])
t.b <- c(t.b, round(tags.a[i]*cons/radios[j]) )
tags <- c(tags, t.a[k]+t.b[k])
r <- c(r, t.a[k]*cons/t.b[k]) # radios REALES
#  tot.col.A
#      To      In 
# 9160694 5743943 
p.A.To <- c(p.A.To, t.a[k]/9160694) 
p.A.In <- c(p.A.In, t.b[k]/5743943)
temp <- forma(a=t.a[k], b=t.b[k]) # Obtenemos la matriz
test <- chisq.test(temp) # Realizamos la prueba
chi <- c(chi, test$statistic) # Capturamos el estad¡stico
P <- c(P, test$p.value) # Capturamos la significancia
}
}

datos <- data.frame(i1, i2, t.a, t.b, tags, r, p.A.To, p.A.In, chi, P)

attach(datos) # Poniendo datos a la vista

# Gr fica
plot(radios, datos$P[datos$i1==1], ylim=c(0,1), type="l", col="black",
xlab = "Radios P[A|To]/P[A|In]", ylab = "Significancia (P)",
main = "Valores de significancia de la prueba chi cuadrada\ncuando se presentan distintos radios y frecuencias de expresi¢n",
sub = "(los radios son aproximados)")
points(radios, datos$P[datos$i1==2], type="l", col="red")
points(radios, datos$P[datos$i1==3], type="l", col="blue")
points(radios, datos$P[datos$i1==4], type="l", col="purple")
points(radios, datos$P[datos$i1==5], type="l", col="orange")
legend("topright", legend=c("1","10","100","1000","10000"),
bty="n", title="P[A|To] en transcritos por millon",
text.col=c("black", "red", "blue", "purple", "orange"))

# prob.A.To
# [1] 1e-06 1e-05 1e-04 1e-03 1e-02
# 1000000*prob.A.To # (expresadas en transcritos por mill¢n)
# [1]     1    10   100  1000 10000

plot(radios[radios<=2], P[(i1==2)&(radios<=2)], ylim=c(0,0.2), type="l", col="red",
xlab = "Radios P[A|To]/P[A|In]", ylab = "Significancia (P)",
main = "Valores de significancia de la prueba chi cuadrada\ncuando se presentan distintos radios y frecuencias de expresi¢n",
sub = "(los radios son aproximados)")
points(radios[radios<=2], P[(i1==3)&(radios<=2)], type="l", col="blue")
points(radios[radios<=2], P[(i1==4)&(radios<=2)], type="l", col="purple")
points(radios[radios<=2], P[(i1==5)&(radios<=2)], type="l", col="orange")
legend("topright", legend=c("10","100","1000","10000"),
bty="n", title="P[A|To] en transcritos por millon",
text.col=c("red", "blue", "purple", "orange"),
title.col="black")

# Acercamiento
plot(radios[radios<=1.2], P[(i1==3)&(radios<=1.2)], ylim=c(0,0.05), type="l", col="blue",
xlab = "Radios P[A|To]/P[A|In]", ylab = "Significancia (P)",
main = "Valores de significancia de la prueba chi cuadrada\ncuando se presentan distintos radios y frecuencias de expresi¢n",
sub = "(los radios son aproximados)")
points(radios[radios<=1.2], P[(i1==4)&(radios<=1.2)], type="l", col="purple")
points(radios[radios<=1.2], P[(i1==5)&(radios<=1.2)], type="l", col="orange")
legend("topright", legend=c("100","1000","10000"),
bty="n", title="P[A|To] en transcritos por millon",
text.col=c("blue", "purple", "orange"),
title.col="black")

# Pruebas m£ltiples
azar <- rnorm(100) # Genero una muestra de 100 n£meros aleatorios de la distribuci¢n normal
summary(azar) # Estad¡sticas b sicas
sd(azar) # Desviaci¢n est ndar
hist(azar) # Obtengo un histograma de los datos

# Muestreando
i1 <- sample(c(1:100), size=50) # 50 n£meros aleatorios muestreados del vector c(1:100)
i2 <- setdiff(c(1:100), i1) # Obtiene el RESTO de los n£meros; i.e., los que no est n en i1
sort(i1); sort(i2) # Para ver como dividimos la muestra

test <- t.test(azar[i1],azar[i2]) # Prueba de t de dos grupos de 50 datos cada uno

test

# Haciendo multiples pruebas...
p.val <- c() # Aqu¡ guardaremos el valor de P.
for( i in 1:1000) {
i1 <- sample(c(1:100), size=50) # 50 n£meros aleatorios muestreados del vector c(1:100)
i2 <- setdiff(c(1:100), i1) # Obtiene el RESTO de los n£meros; i.e., los que no est n en i1
test <- t.test(azar[i1], azar[i2])
p.val <- c(p.val, test$p.value) # Capturamos el valor de P
}

summary(p.val)
length(p.val[p.val<0.05]) # N£mero de elementos de p.val menores de 0.05 (debe ser m s o menos 50, es decir el 0.05 (5%) de los casos).

length(p.val[p.val<5e-05]) # N£mero de elementos de p.val menores de 5e-05



9