RESUM de funcions amb R (EXAMEN FINAL) (2014)

Resumen Catalán
Universidad Universidad de Barcelona (UB)
Grado Biotecnología - 4º curso
Asignatura Disseny Experimental i Anàlisi de Dades (DEAD)
Año del apunte 2014
Páginas 9
Fecha de subida 30/03/2016
Descargas 14
Subido por

Vista previa del texto

Més apunts a Unybook.com, buscant l’usuari: biotechub options(show.signif.stars=FALSE, contrasts=c("contr.sum", "contr.poly")) t.test(x,y,var.equal=TRUE) - 'y' s'omet si es tracta només una mostra (llavors la variable d'estudi ja és la diferència de mitjes en si).
- Modificar mu (mu=el que sigui).
- Modificar H1 (alternative="greater" o bé alternative="less") - Dades aparellades (mostres dependents) (paired=TRUE).
qt(0.975,df=n+m-2) qn qf (0.975,df1=4,df2=20) especificar els dos graus de llibertat var.test(x,y) chisq.test(x) on 'x' serà una taula (de contingència) Test t de Student de comparació de mitjanes de 'x' i 'y'. var.equal=TRUE indica que assumim que les variàncies són iguals. Entre d'altres, retorna el valor de l'estadístic 't', el p-valor i l'interval de confiança del 95%. Per defecte considera una hipòtesi alternativa bilateral (µ1≠µ2). Per veure tots els paràmetres que vols incloure, prémer tabulador dins del parèntesi (això serveix per totes les funcions). Si es tracta de mostres dependents (ex: colesterol abans i després de una dieta) => comparació de mitjanes per mostres normals aparellades (paired=TRUE), en aquest cas no cal test de comparació de variàncies ja que són iguals.
Troba el valor crític d'una t-Student amb n+m-2 graus de llibertat (comparació de mitjanes de dues poblacions assumint igualtat de variàncies). 0'975 => nivell de significació de 0'05 i hipotesi alternativa bilateral.
El mateix que a dalt però per una distribució Normal El mateix que a dalt però per una distribució F de Fisher Test t F de Fisher de comparació de variàncies de 'x' i 'y'. Entre d'altres, retorna el valor de l'estadístic 't', el p-valor i l'interval de confiança del 95%.
Test khi quadrat. S'utilitza per tests d'homogeneïtat i d'independència.
Especificar 'correct=TRUE' si apliques la correcció per continuïtat (s'aplica quan alguna de les freqüències de la taula és inferior a 5). Les dades han d'estar en forma de taula, per fer-ho: 1) Introduir el vector de dades, 2) Definir l'objecte com a matriu aplicant matrix(objecte,nrow=A,ncol=B) on A i B son el nº de files i columnes respectivament, i 3) Definir la matriu com a taula aplicant as.table(objecte). Apliques la Khi Quadrat al nou objecte que ha quedat redefinit com a taula.
Llegeix un .txt i si posem header=TRUE vol dir que en el .txt has posat una capçalera/títol a les columnes.
read.table("doc.txt", header=T) ex: porphyra <- read.table("PORPHYRA.DAT", header = TRUE) Datos$grupo <Etiquetar nivells. Això ho fem pq R agafi grupo com una variable factor, sinó ho agafa com factor(Datos$grupo,labels=c('15%','20%','25%','30%','35%')) una variable numèrica. Les cometes indiquen que és text! Com a primer argument passem una fórmula. En R un objecte de classe 'formula' sempre té la forma 'variable de resposta' ~ 'expressió'. A 'expressió' hi intervenen diferents variables, a Més apunts a Unybook.com, buscant l’usuari: biotechub anova(lm(medida ~ grupo, data=Datos)) , on medida és la variable resposta i grupo és factor Coto_a.aov <- aov(medida ~ grupo, data = Datos) summary(Coto_a.aov) model.tables(Coto_a.aov,type="mean") tapply(Datos$medida, Datos$grupo, mean, na.rm=TRUE) tapply(Datos$medida, Datos$grupo, sd, na.rm=TRUE) tapply(Datos$medida, Datos$grupo, function(x) sum(!is.na(x))) tapply(Datos$medida, Datos$grupo, median, na.rm=TRUE) library(Rcmdr) plotMeans (Datos$medida, Datos$grupo, error.bars="conf.int", level=0.95) intervals confiança 95% dummy.coef(Coto_a.aov) residuals(Coto_a.aov) VALIDESA DEL MODEL: els errors han de verificar les condicions de gauss-markov plot(Coto_a.aov) qqnorm(residuals(Coto_a.aov)) boxplot(medida~grupo, ylab="medida", xlab="grupo", data=Datos) boxplot(Esporangis ~ Color, data = PORPHYRA) partir de les quals voldríem explicar els valors de la variable de resposta. Com a segon argument indiquem a quin 'data.frame' hi ha les variables que corresponen als noms de variable que s'han fet servir a la fórmula.
ANOVA Mostra la taula ANOVA, conté MSe (l'arrel d'això és la sigma estimada) Taula de mitjanes: (mu's amb sombrero): mitjanes estimades de cada nivell Mitjana Desviació estandard Tamany dels nivells Mediana tapply: permet aplicar una funció per a cada nivell del factor.
El paquet Rcmdr permet fer gràfics com: representació de les mitjanes (plotMeans) Coeficients del model (estimació dels efectes alfa de cada nivell), intercept = estimació mu general Estimació dels errors VALIDESA DEL MODEL: els errors han de verificar les condicions de gauss-markov Si els residus apareixen raonablement alineats => suggereix que segueix normalitat Veure observacions (respostes) extremes.
Homocedasticitat: homogeneitat de variança Barlett no és molt bo si no es compleix normalitat.
Levene es mes robust si no es compleix normalitat (si es desvia). Levene requereix la llibreria lawstat.
Més apunts a Unybook.com, buscant l’usuari: biotechub tapply(Datos$medida, Datos$grupo, var, na.rm=TRUE) bartlett.test(medida ~ grupo, data=Datos) resposta ~factor, data= x require(lawstat) levene.test(Datos$medida, Datos$grupo) x$resposta, x$factor kruskal.test(medida ~ grupo, data=Datos) TEST DE COMPARACIONS MÚLTIPLES (test de mitjanes a posteriori) (si l'ANOVA diu q hi ha diferencies significatives) require(agricolae) LSD.test(Coto_a.aov, "grupo", p.adj=c("none")) LSD.test(Coto_a.aov, "grupo", p.adj=c("bonferroni")) No es suposa normalitat ni homocedasticitat. ALTERNATIVA NO PARAMÈTRICA AL ANOVA Test LSD (Least Significant Difference = valor crític del test) LSD.test(Coto_a.aov, "grupo", p.adj=c("none"), group=F) LSD.test(Coto_a.aov, "grupo", p.adj=c("bonferroni"), group=F) Per obtenir les probabilitats de les comparacions s'ha de indicar que no es requereix grups.
# Com que LSD realment no controla l'error de tipus I, cal especificar un # mètode addicional de correcció per multiplicitat (opcions: "none","holm", # "hochberg", "bonferroni", "BH", "BY", "fdr") # Una bona opció podria ser "holm" que controla l'error de tipus I global i # és vàlid per p-valors dependents entre ells.
LSD.test(porphyra.aov, "Color", p.adj = "holm") LSD.test(porphyra.aov, "Color", p.adj = "holm", group = FALSE) duncan.test(Coto_a.aov, "grupo", group=F) Prova de Duncan. Nivell de significació real més gran que alfa (al risc error tipus I) però també alta potència.
SNK.test(Coto_a.aov, "grupo", group=F) Prova SNK (Student-Newman-Keuls) Més apunts a Unybook.com, buscant l’usuari: biotechub HSD.test(Coto_a.aov, "grupo", group=F) TukeyHSD(porphyra.aov) (aquesta esta en el paquet base) Prova de TukeyHSD Diff = diferencia de mitjanes. Dóna directament el p-valor ajustat O també: require(multcomp) porphyra.glht <- glht(porphyra.aov, linfct = mcp(Color = "Tukey")) summary(porphyra.glht) confint(porphyra.glht) plot(porphyra.glht) El plot mostra un interval de confiança del 95%. Aquest plot mostra les diferències entre els nivells del factor (les mes extremes són diferències mes grans) (si el interval inclou el zero NO hi ha diferències significatives) require(multcomp) porphyra.glht <- glht(porphyra.aov, linfct = mcp(Color = "Dunnett")) summary(porphyra.glht) confint(porphyra.glht) plot(porphyra.glht) Prova de Dunnet Per definir un nivell específic com a control: (en aquest cas: 'roja' porphyra.glht <- glht(porphyra.aov, linfct = mcp(Color = c("roja - blava = 0", "roja - verda = 0", "roja - groga = 0", "roja - infra = 0"))) on Color és el factor, i els diferents colors són nivells.
La funció confint() mostra els intervals confiança 95% per la diferència de mitjanes.
scheffe.test(porphyra.aov, "Color") scheffe.test(porphyra.aov, "Color", group = FALSE) O bé require(multcomp) porphyra.glht <- glht(porphyra.aov, linfct = mcp(Color = "Schefe")) summary(porphyra.glht) confint(porphyra.glht) plot(porphyra.glht) Prova de Schéffe Més apunts a Unybook.com, buscant l’usuari: biotechub DISSENY DE DOS FACTORS AMB INTERACCIÓ # Els factors convé que R els entengui com variables de classe "factor" (no numèrica, per Per saber quin tipus de dada és el que hem entrat s'utilitza la exemple). Aquí no cal, però si algun dels factors anteriors no els considerés com a tals, funció class() caldria fer: blat$Fertilitzant <- as.factor(blat$Fertilitzant).
blat.aov <- aov(Produccio ~ Fertilitzant * Varietat, data = blat) summary(blat.aov) per obtenir la taula anova.
dummy.coef(blat.aov) residuals(blat.aov) Fer gràfic QQ plot(residuals(blat.aov)) per veure si els errors segueixen una normal i si hi ha punts atípics (que també es poden detectar amb boxplot). També es pot fer qqnorm. Per fer el boxplot dels residus: e <- residuals(blat.aov) boxplot(e) Fer l'ANOVA. # En un model lineal especificat en R, el símbol '*' vol dir "el model creuat # complet", és a dir, el model que inclou l'efecte per separat de cada un # dels factors i la interacció.
# Seria equivalent a: #La interaccio s'escriu amb el nom de cada factor separat per dos punts (:).
blat.aov <- aov(Produccio ~ Fertilitzant + Varietat + Fertilitzant:Varietat, data = blat) # El símbol ':' permet especificar una interacció. Així amb l'ordre: blat.aov <- aov(Produccio ~ Fertilitzant + Varietat, data = blat) # ajustaríem un model sense interacció. O amb blat.aov <- aov(Produccio ~ Fertilitzant + Fertilitzant:Varietat, data = blat) # ajustaríem un model sense l'efecte separat de 'Varietat', però amb interacció.
Estimacions del coeficients del model lineal Produccio ~ Fertilitzant * Varietat.
Estimacions dels residus (errors).
# Aquesta funció dóna una informació similar que dummy.coef: model.tables(blat.aov) #proporciona la estimacio dels alfes, betes, gammes, de manera organitzada.
# Si volem que per cada estimació ens mostri també el seu error estàndard: model.tables(blat.aov, se=TRUE) # Si en lloc dels coeficients del model volem les taules de mitjanes: # separadament per cada nivell dels factors i per cada combinació de nivells o # condició experimental ("casella"): model.tables(blat.aov, type="means", se=TRUE) # A la pregunta: "indica totes les estimacions dels paràmetres del model", # la resposta més correcta seria indicar els valors retornats per 'dummy.coef' # (o per 'model.tables'), i el valor de la variancia residual que podem treure de summary.
Més apunts a Unybook.com, buscant l’usuari: biotechub interaction.plot(blat$Varietat, blat$Fertilitzant, blat$Produccio) #L'ordre és important (variable resposta va l'últim). Amb la interaction plot ja podem escollir quina és la millor combinacio.
interaction.plot(blat$Fertilitzant, blat$Varietat, blat$Produccio, col=1:length(levels(blat$Varietat))) # Per representar 'Fertilitzant' a l'eix d'abcisses i donar un color diferent # a cada línia (cada línia representa un nivell de 'Varietat') model.tables(blat.aov, type="means", cterms = "Fertilitzant:Varietat", se=TRUE) # Observació de les mitjanes per cel·les (o "caselles" o "condicions # experimentals, és a dir, combinacions dels nivells de cada factor), que és # la que té més sentit sota interacció. Per trobar la condició experimental òptima.
COMPARACIONS MÚLTIPLES TukeyHSD(blat.aov, which="Fertilitzant") comp.mult <- TukeyHSD(blat.aov, which="Fertilitzant") plot(comp.mult) blatA.aov <- aov(Produccio ~ Fertilitzant, data = blat, subset = blat$Varietat == "A") # == implica comparacio.
comp.multA <- TukeyHSD(blatA.aov, which="Fertilitzant") #Dins de la varietat A quines diferencies hi ha entre fertilitzants comp.multA #Com es veu, dins la varietat A tots són diferents plot(comp.multA) inter.aov = aov(Produccio ~ Fertilitzant:Varietat, data = blat) #és un anova d'un factor, però el factor es la interaccio model.tables(inter.aov, type="means", cterms = "Fertilitzant:Varietat", se=TRUE) comp.multI <- TukeyHSD(inter.aov) comp.multI Quan hi ha interacció, fer comparacions múltiples marginalment per cada factor no té gaire sentit. Si de tota manera les volguéssim fer pel factor Fertilitzant (ja que hem vist que hi ha diferències significatives per algun d'ells).
Per fer un plot amb l'interval de confiança de la diferència de mitjanes dels fertilitzants.
# Però tal com hem dit, en aquest cas que hi ha evidències d'interacció, no és # recomanable fer comparacions múltiples de "Fertilitzant", per molt que hagi # sortit significatiu.
# Normalment, en aquest cas faríem comparacions entre fertilitzants PER SEPARAT # DINS CADA NIVELL DE L'ALTRE FACTOR (o dins els nivells de "Varietat" que ens # semblessin més interessants).
# Encara que possiblement la comparació més interessant, atès que hi ha # interacció significativa, és entre les cel·les (combinacions dels nivells # de cada factor): Més apunts a Unybook.com, buscant l’usuari: biotechub homogeneïtat de variàncies: library(car) # Entre condicions experimentals: leveneTest(Produccio ~ Fertilitzant*Varietat, data = blat) # Entre nivells de Fertilitzant: leveneTest(Produccio ~ Fertilitzant, data = blat) FUNCIÓ attach(document amb les dades) permet estalviar-te haver d'escriure els $.
Test de homogeneitat -> test Khi quadrat per taules de contingència. Per introduir les dades d'una taula -> fer-ho per columnes. Ex: freq <- c(100,71,29,146,32,17,149,28,16,146,37,17,147,33,18) Ara hem de transformar el vector en matriu, fent servir matrix() i indicant el numero de files amb 'nrow' (el numero de columnes es dedueix a partir del de files). Tambe podriem posar ncol=5 freq <- matrix(freq, nrow=3) El seguent pas amb R seria convertir aquesta matriu en una taula. No es necessari per fer el test, però si per ser coherent (R es un llenguatge de programacio orientat a objectes, és a dir, depenent del tipus de objectes fa unes coses o unes altres. P.ex: una table té associada certes funcions, mentre que una matric en té associades unes altres) freq <- as.table(freq) Ara R sap que es una taula de freq. i podem aplicar-hi el test de la Khi quadrat.
Més apunts a Unybook.com, buscant l’usuari: biotechub Introducció de dades #Es fa un vector amb totes les mesures i un altre que et classifica les dades. rep = repetir. a=actual, f=fossil longitud <- c(0.92,1.29,1.0,1.5,1.25,1.54,1.26,1.71,1.28,1.16,1.43,1.8,1.5,1.57,1.75) poblacio <- c(rep("a",9),rep("f",6)) #és el mateix que: rep(c("a","f"),c(9,6)) #Hem de convertir el caràcter poblacio a factor. La funcio data.frame ho fa directament: transforma vector caràcter a vector factor.
siurana <- data.frame(longitud,poblacio) #Ara el data.frame conté els dos vectors. Per tant tenim tres objectes (siurana, longitud i població).
#En la següent instruccio (rm=remove) eliminem els vectors longitud i poblacio perque R no s'equivoqui i em calculi les funcions en un objecte que no m'interessa.
rm(longitud,poblacio) siurana summary(siurana) #attach es una instrucció de R que ens permetrà utilitzar els noms poblacio i longitud sense el dolar (el dolar ens serveix per indicar on esta el vector, ex: poblacio dintre de siurana -> poblacio$siurana). Si abans de fer attach posem 'longitud' a la consola ens diu que no el troba (perquè està dintre de siurana) attach(siurana) #al fer attach i tornar a escriure 'longitud' a la consola veiem com si que el troba tapply(longitud, poblacio, mean) #la 't' de tapply ve de 'tractament'. Significat d'aquest 'tapply': al vector longitud li calculem la mitjana separat per poblacio. Una altra manera de calcular la mitjana per separat (per cada poblacio): mean(longitud[poblacio=="a"]) mean(longitud[poblacio=="f"]) Més apunts a Unybook.com, buscant l’usuari: biotechub fercoefs <- rep(c(0,21.5,50.5), 4) #es repeteix 4 vegades pq hi ha 4 varietats. Executa rep i obtenim les estimacions varcoefs <- rep(c(0,5,-5.3,6),each=3) #each = repetir seguit. Executa rep i obtenim les estimacions Pas de vector numèric a vector factor: pes_nado <- c(3.205, 4.10, 3.55, 2.30, 4.30,3.435,3.30,2.560, 3.600,4.560,3.440,4.340,3.560,4.550,2.550) # Creació d' un objeto 'Tipo Mare' que codifica la mare para cada uno de los valores anteriores: tipo_mare <- factor(rep(c(1,2,3,4,5), c(3,3,3,3,3))) --Iris.setosa <- c(5.1, 4.9, 4.7, 4.6, 5, 5.4, 4.6, 5, 4.4, 4.9, 5.4, 4.8, 4.8, 4.3, 5.8) Iris.versicolor <- c(7, 6.4, 6.9, 5.5, 6.5, 5.7, 6.3, 4.9, 6.6, 5.2, 5, 5.9, 6, 6.1, 5.6) Iris.virginica <- c(6.3, 5.8, 7.1, 6.3, 6.5, 7.6, 4.9, 7.3, 6.7, 7.2, 6.5, 6.4, 6.8, 5.7, 5.8) #Pero esta no es la forma adecuada para trabajar con un programa estadistico. Mejor ponemos los datos en un unico vector y a~nadimos una variable cualitativa o factor que nos indique la poblacion de cada dato.
longitud <- c(Iris.setosa, Iris.versicolor, Iris.virginica) #Fent-ho aixi puc combinar els tres vectors en un de sol especie <- rep(1:3, each = 15) #és el mateix que posar: 1,1,1,1,1 (15 vegades), 2,2,2,2,2 (15 vegades), 3,3,3,3,3 (15x) #O també podriem escriure: # especie <- c(rep(1,15),rep(2,15), rep(3,15)) #Amb qualsevol de les maneres aconseguirem repartir les dades en 3 columnes. Ara passem especie a variable factor especie <- factor(especie, labels = c("Iris setosa", "Iris versicolor","Iris virginica")) ...