Matemàtiques II Graus de Biologia i Bioquímica 13-14 UIB
La taula de dades algae de R consisteix en mostres d’aigua preses en distints rius europeus agafades en distintes estacions de l’any. Cada observació (mostra d’aigua) té variables. Tres d’aquestes variables són nominals i descriuen l’estació de l’any quan les mostres d’aigua varen ésser preses (season), la grandària (size) i la velocitat del riu (speed). Les vuit variables següents són mesures de paràmetres químics:
Les variables restants són el nombre d’algues nocives que es troben a la mostra d’aigua respectiva. N’hi ha perquè s’han considerat tipus d’algues nocives. Per carregar les dades heu d’anar a http://bioinfo.uib.es/~recerca/MAT2/tallers/algae.txt
Escriviu una funció que us doni l’interval de confiança per a la proporció per una mostra de grandària d’una variable de Bernoulli amb proporció d’èxit desconeguda. Els paràmetres d’entrada de la funció han d’ésser (mostra) i (nivell de significació). La sortida ha d’ésser un vector amb dues components: els dos extrems de l’interval de confiança per a la proporció. La vostra funció ha de tenir en compte que la mostra pot tenir dades perdudes (NA) en el sentit que les ha d’eliminar. Poseu també que el valor per defecte de sigui La vostra funció també ha de tenir en compte que si la mostra no és una mostra corresponent a una variable de Bernoulli, ha d’avisar i escriure el missatge següent: “La mostra no correspon a una variable de Bernoulli”. Digueu-li a la funció IntConfProporcio.
Recordem que si \( x_1,x_2,...,x_n \) es una realització d'una m.a.s. \( Ber(p) \) llavors
\( \hat{p}=\frac{\sum_{i=1}^n x_i}{n} \) es la estimació de la proporció mostral i el interval de confiança al nivell \( (1-\alpha)\cdot 100\% \) es pot aproximar per mostres grans per
\[ \left(\hat{p}+z_{\alpha/2}\cdot \sqrt{\frac{\hat{p}\cdot (1-\hat{p})}{n}},\hat{p}+z_{1-\alpha/2}\cdot \sqrt{\frac{\hat{p}\cdot (1-\hat{p})}{n}}\right). \]
IC.Ber = function(mostra, alpha) {
mostra = na.omit(mostra)
n = length(mostra)
prop.muestral = sum(mostra)/n
IC = c(prop.muestral + qnorm(alpha/2) * sqrt(prop.muestral * (1 - prop.muestral)/n),
prop.muestral + qnorm(1 - alpha/2) * sqrt(prop.muestral * (1 - prop.muestral)/n))
return(IC)
}
# Exemple
mostra = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0)
Interval0.05 = IC.Ber(mostra, alpha = 0.05)
Interval0.05
## [1] 0.2964 0.9036
Considerem la variable que ens dóna el valor mitjà de clorurs de cada mostra. Construïu quatre variables de Bernoulli de la forma següent: pren el valor si el valor de clorurs per la mostra corresponent agafada en primavera és superior (estricte) a 40 i en cas contrari, el mateix per les mostres agafades a l’estiu, tardor i hivern. Anomeneu ClBPrimavera, ClBEstiu, ClBTardor, ClBHivern a aquestes variables.
La variable dels clorurs te per etiqueta Cl
# Carregam algae
algae = read.table("http://bioinfo.uib.es/~recerca/MAT2/tallers/algae.txt",
header = TRUE, sep = " ")
ClB = algae$Cl
ClB = as.integer(ClB > 40)
ClB
## [1] 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0
## [24] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1
## [47] 1 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA 0 0 0 0 0 0 0
## [70] 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 1
## [93] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1
## [116] 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0
## [139] 0 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 NA
## [162] 1 1 1 0 0 1 1 1 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
## [185] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 NA 1
ClBPrimavera = ClB[algae$season == "spring"]
ClBPrimavera
## [1] 1 1 0 0 1 0 0 0 0 0 1 1 0 0 NA NA NA 0 0 0 1 0 0
## [24] 0 1 1 1 1 0 1 1 1 0 0 0 1 0 0 1 1 1 1 1 NA 1 1
## [47] 1 1 0 0 0 0 0
ClBEstiu = ClB[algae$season == "summer"]
ClBEstiu
## [1] 1 0 0 0 0 0 0 1 0 NA NA 0 0 0 0 0 0 1 1 1 1 1 1
## [24] 0 0 1 0 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1
ClBTardor = ClB[algae$season == "autumn"]
ClBTardor
## [1] 1 1 1 0 0 0 0 0 1 0 NA NA 0 0 0 0 0 1 1 1 0 1 0
## [24] 0 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 1
ClBHivern = ClB[algae$season == "winter"]
ClBHivern
## [1] 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 0 0 0 NA 0
## [24] 0 1 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 1 0
## [47] 1 1 1 1 0 0 1 1 1 1 0 0 0 0 0 NA
Continuem amb la variable de l’exercici 2. Calculeu ara els intervals de confiança per a la proporció de mostres a un 95% de confiança on el seu nivell de clorurs sigui superior a 40 per mostres agafades a la primavera, estiu, tardor i hivern. Hi ha qualque indicació de que en qualque estació hi hagi més mostres amb un nivell de clorur superior a 40?
IC.Ber(ClBPrimavera, 0.05)
## [1] 0.3297 0.6091
IC.Ber(ClBEstiu, 0.05)
## [1] 0.2492 0.5415
IC.Ber(ClBTardor, 0.05)
## [1] 0.2641 0.5780
IC.Ber(ClBHivern, 0.05)
## [1] 0.3241 0.5759
La primavera es la estació amb el interval més alt.
Considerem ara la variable de Bernoulli ClB que ens dóna 1 si el nivell de clorur d’una mostra determinada és superior (estricte) a 40 i 0 en cas contrari. Calculeu 100 mostres de tamany 100 amb reposició de la variable anterior. Guardau aquestes mostres en una matriu anomenada mostres de 100 files i 100 columnes. Pensau a llevar prèviament tots els valors NA de la variable clorur Cl.
ClB = na.omit(ClB)
mostra = matrix(sample(ClB, size = 100^2, replace = TRUE), nrow = 100, ncol = 100)
class(mostra)
## [1] "matrix"
dim(mostra)
## [1] 100 100
head(mostra, 3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 0 1 0 0 1 0 1 1 1 1 1 1 1
## [2,] 0 1 1 1 1 1 1 1 0 1 0 1 0
## [3,] 1 0 0 0 0 1 0 0 0 0 1 1 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 0 0 1 0 1 0 1 0 0 1 1
## [2,] 1 0 1 0 0 0 1 0 1 0 1
## [3,] 1 0 1 0 0 1 1 1 0 0 0
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 1 1 0 1 0 1 0 0 0 1 1
## [2,] 1 0 0 0 1 1 1 1 1 1 0
## [3,] 0 0 1 1 0 0 1 0 0 0 1
## [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,] 0 1 1 0 0 1 0 1 0 0 0
## [2,] 1 0 1 0 1 0 1 0 0 1 0
## [3,] 1 1 0 1 0 0 1 1 0 1 1
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57]
## [1,] 0 0 1 1 0 1 0 1 1 0 0
## [2,] 1 0 0 1 1 0 1 0 0 0 0
## [3,] 1 1 1 1 0 1 0 0 0 0 0
## [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,] 0 0 0 1 1 0 1 0 0 0 1
## [2,] 0 0 0 0 0 1 0 1 0 0 0
## [3,] 0 0 1 0 1 0 0 0 1 1 0
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79]
## [1,] 0 0 0 1 1 0 0 1 0 0 0
## [2,] 0 1 0 1 1 0 0 1 1 0 1
## [3,] 0 0 1 0 0 1 1 0 0 1 0
## [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] 0 0 1 0 0 0 1 0 0 0 0
## [2,] 0 0 0 0 0 1 0 0 0 1 1
## [3,] 1 1 0 0 0 0 0 1 1 0 0
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## [1,] 0 0 1 0 1 1 1 0 1 0
## [2,] 1 1 0 0 1 0 0 1 0 1
## [3,] 0 0 1 1 1 1 1 0 1 1
Fent servir la funció definida en l’exercici 1, calculeu els corresponents intervals de confiança a un nivell del de confiança per la proporció i compteu en quants d’aquests intervals estarà la proporció de mostres de la variable clorur superior a 40. Quin és el nombre esperat d’intervals en els que la proporció de mostres de la variable clorur és superior a 40 hi estarà dins? Coincideixen els valors obtinguts amb els valors esperats? Raoneu les respostes.
proporcio = sum(ClB)/length(ClB)
proporcio
## [1] 0.4368
IC100Mostres = apply(mostra, MARGIN = 1, FUN = IC.Ber, alpha = 0.05)
IC100Mostres
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.333 0.3623 0.3427 0.3623 0.304 0.3525 0.3722 0.4221 0.3623 0.3233
## [2,] 0.527 0.5577 0.5373 0.5577 0.496 0.5475 0.5678 0.6179 0.5577 0.5167
## [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 0.304 0.2659 0.2944 0.392 0.304 0.3136 0.2849 0.3136 0.392 0.3525
## [2,] 0.496 0.4541 0.4856 0.588 0.496 0.5064 0.4751 0.5064 0.588 0.5475
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
## [1,] 0.3525 0.2944 0.392 0.304 0.3623 0.2565 0.333 0.3623 0.304 0.3722
## [2,] 0.5475 0.4856 0.588 0.496 0.5577 0.4435 0.527 0.5577 0.496 0.5678
## [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40]
## [1,] 0.3525 0.333 0.3525 0.3821 0.402 0.2849 0.2472 0.3623 0.2659 0.3233
## [2,] 0.5475 0.527 0.5475 0.5779 0.598 0.4751 0.4328 0.5577 0.4541 0.5167
## [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 0.402 0.392 0.304 0.473 0.3136 0.3233 0.3623 0.3722 0.3525 0.2944
## [2,] 0.598 0.588 0.496 0.667 0.5064 0.5167 0.5577 0.5678 0.5475 0.4856
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
## [1,] 0.3525 0.2849 0.392 0.3427 0.333 0.2472 0.333 0.3623 0.3623 0.2849
## [2,] 0.5475 0.4751 0.588 0.5373 0.527 0.4328 0.527 0.5577 0.5577 0.4751
## [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70]
## [1,] 0.333 0.2565 0.3821 0.3136 0.3136 0.2659 0.3427 0.2944 0.3427 0.3233
## [2,] 0.527 0.4435 0.5779 0.5064 0.5064 0.4541 0.5373 0.4856 0.5373 0.5167
## [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [1,] 0.3136 0.304 0.3427 0.4221 0.333 0.2944 0.2659 0.3623 0.304 0.3821
## [2,] 0.5064 0.496 0.5373 0.6179 0.527 0.4856 0.4541 0.5577 0.496 0.5779
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] 0.3821 0.392 0.304 0.304 0.3136 0.3623 0.3427 0.2754 0.3722 0.3136
## [2,] 0.5779 0.588 0.496 0.496 0.5064 0.5577 0.5373 0.4646 0.5678 0.5064
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## [1,] 0.3136 0.2849 0.3233 0.333 0.412 0.333 0.3623 0.412 0.2659 0.2754
## [2,] 0.5064 0.4751 0.5167 0.527 0.608 0.527 0.5577 0.608 0.4541 0.4646
aux = IC100Mostres[1, ] < proporcio & IC100Mostres[2, ] > proporcio
aux
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [12] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [23] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [34] TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [45] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [56] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [67] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [78] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [89] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [100] TRUE
table(aux)
## aux
## FALSE TRUE
## 3 97
Un 97 % dels intervals de confiança contenen la proporció mostral. Teòricament el 95% dels intervals contendrien la proporció muestral.
El resultat es pot observar en el següent gràfic
require(Hmisc)
## Loading required package: Hmisc Loading required package: survival Loading
## required package: splines Loading required package: Formula Hmisc library
## by Frank E Harrell Jr
##
## Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview') to see overall
## documentation.
##
## Attaching package: 'Hmisc'
##
## The following object is masked from 'package:survival':
##
## untangle.specials
##
## The following object is masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, units
Mostres = c(1:100)
xYplot(Cbind((IC100Mostres[1, ] + IC100Mostres[2, ])/2, IC100Mostres[1, ], IC100Mostres[2,
]) ~ Mostres, groups = c(1:100), hotizontal = TRUE, col = c("red", "blue")[unclass(as.factor(aux))],
ylab = "Interval de confiança al 95%")
## Loading required package: grid Loading required package: lattice