#source("/Users/Emmanuel/Documents/Outils/R-Fonctions.R") LastChar <- function (x) { L = nchar(x) return(substr(x,L,L)) } Detacher <- function (nom) { S <- search() while ( length(grep(nom, S))>0 ) { detach(pos=grep(nom, S) [1]) S <- search() } } meanNA <- function (x) {mean(x,na.rm=T)} ebars <- vbars <- function (x, y0, y1, y2, ...) { for(i in 1:length(x)) { if(y0[i]!=0 & !is.na(y0[i]) & y1[i]!=0 & !is.na(y1[i]) & (y0[i]!=y2[i])) { arrows (x[i], y0[i], x[i], y1[i], angle=90, length=0.05, ...) } if(y0[i]!=0 & !is.na(y0[i]) & y2[i]!=0 & !is.na(y2[i]) & (y0[i]!=y2[i])) { arrows (x[i], y0[i], x[i], y2[i], angle=90, length=0.05, ...) } } } StdErr <- function(x) { try(return(sqrt(var(x,na.rm=TRUE)/length(x))), silent=T) return(10) } BarErr1 <- function(data, niveau, restrictions = TRUE, ...) { xM <- tapply(data[restrictions], niveau[restrictions], mean, na.rm=T) xS <- tapply(data[restrictions], niveau[restrictions], StdErr) x <- barplot(xM, ...) ebars(x, xM, xM+xS, xM-xS) } BarErr2 <- function(data, niveau1, niveau2, restrictions = TRUE, ...) { xM <- tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), mean, na.rm=T) xS <- tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), StdErr) x <- barplot(xM, ...) ebars(x, xM, xM+xS, xM-xS) } LigneLateX <- function(A) { cond=A[1] val1=A[2] dev1=A[3] val2=A[4] dev2=A[5] return(c("\\drawcondB","{",cond,"}{",val1,"}{",dev1,"}{",val2,"}{",dev2,"}\\\\\n", sep="")) } GraphLateX <- function(data, niveau1, niveau2, colonnes, restrictions = TRUE, ...) { xM <- round(tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), mean, na.rm=T),2) xS <- tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), StdErr) d <- merge(xM, xS, by="row.names") e <- c("Row.names",paste(rep(colonnes,each=2), c(".x", ".y"), sep="")) RES <- c("\\fbox{\\begin{tabular}{r|rl|rl}\n & \\multicolumn{2}{c|}{", colonnes[1], "}", "& \\multicolumn{2}{c}{", colonnes[2], "}\\\\\n") RES <- c(RES, apply(d[e], 1, LigneLateX) ) RES <- c(RES, "\\end{tabular}}") return(RES) } LineErr2 <- function(data, niveau1, niveau2, restrictions=TRUE, ...) { xM <- tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), mean, na.rm=T) xS <- tapply(data[restrictions], list(niveau1[restrictions], niveau2[restrictions]), StdErr) # x <- barplot(xM, col="white", ...) plot(x, type="n", ylim=range(0, xM+xS, xM-xS), xlim=range(1,length(xM[1,])), ...) lines(1:length(xM[1,]), xM[1,], col="blue") lines(1:length(xM[2,]), xM[2,], col="red") ebars(1:length(xM[1,]), xM[1,], xM[1,]+xS[1,], xM[1,]-xS[1,], col="blue", lwd=2) ebars(1:length(xM[2,]), xM[2,], xM[2,]+xS[2,], xM[2,]-xS[2,], col="red", lwd=2) points(matrix(rep(1:length(xM[1,]),2),nrow=2, byrow=T), xM, col=c("blue", "red"), pch=12) } # Plus simple: # N <- 40 # A <- rnorm(N) # N1 <- 5 # N2 <- 2 # F1 <- gl(N1,1,N) # F2 <- gl(N2,3,N) LineErr <- function(data, niveau1, niveau2, restriction = TRUE, ...) { N1 <- length(levels(niveau1[restriction])) N2 <- length(levels(niveau2[restriction])) N <- length(data[restriction]) dessin=1:N2 M <- matrix(rep(1:N1,N2), nrow=N2, byrow=T) xM <- tapply(data[restriction],list(niveau2[restriction], niveau1[restriction]),mean, na.rm=T) xS <- tapply(data[restriction],list(niveau2[restriction], niveau1[restriction]),StdErr) interaction.plot(niveau1[restriction], niveau2[restriction], data[restriction], fun=meanNA, type="b", lty=dessin, col=dessin, pch=dessin, ...) arrows(M, xM-xS, M, xM+xS, code=3, angle=90, length=.05, col=dessin)#, ...) #MARCHE PAS BIEN: #legend("topleft", levels(F2), pch=1:N2+15, col=1:N2, lty=1:N2, cex=2) } LineErr <- function(data, niveau1, niveau2, restriction = TRUE, ...) { N1 <- length(levels(niveau1[restriction])) N2 <- length(levels(niveau2[restriction])) N <- length(data[restriction]) dessin=1:N2 M <- matrix(rep(1:N1,N2), nrow=N2, byrow=T) xM <- tapply(data[restriction],list(niveau2[restriction], niveau1[restriction]),mean, na.rm=T) xS <- tapply(data[restriction],list(niveau2[restriction], niveau1[restriction]),StdErr) interaction.plot(niveau1[restriction], niveau2[restriction], data[restriction], fun=meanNA, type="l", lty=1, col=dessin, pch=0, ...) arrows(M, xM-xS, M, xM+xS, code=3, angle=90, length=.05, col=dessin)#, ...) #MARCHE PAS BIEN: #legend("topleft", levels(F2), pch=1:N2+15, col=1:N2, lty=1:N2, cex=2) } EnleverSD <- function(x, y, r) { #renvoie les "indices" des donnees de x a moins de 2 ecarts #types within y[r] #r est une preselection des donnees interessantes (ex: item) xr <- x[r] yr <- y[r] Ecart <- tapply(xr, yr, sd, na.rm=T) Moy <- tapply(xr, yr, mean, na.rm=T) res <- c() for (i in 1:length(x)) { if (r[i] && abs(x[i]-Moy[[y[i]]])<2*Ecart[[y[i]]]) { res <- c(res, TRUE) } else { res <- c(res, FALSE) } } res } NAdonneval <- function(x, val) { if (is.na(x)) {val} else {x} } EnleverNA <- function(x) { #renvoie x en remplacant tous les NA par la moyenne de x apply(x, MARGIN=c(1,2), FUN=NAdonneval, val=mean(x, na.rm=TRUE)) } library(sfsmisc) ###### From Vincent de Gardelle # pour sortir les resultats de l'anova au format standard myprintaov <- function(a, maineff="off") { # a est un objet anova, ex: a <- aov(rt~condition1 + Error(suj/condition1),data=priming) err_lay <- names(summary(a)) for (ierr in err_lay) { cur_lay <- eval(parse(text=paste("summary(a)$ \"",ierr,"\"",sep="")))[[1]] nb_fact_lay <- nrow(cur_lay)-1 nm_fact_lay <- rownames(cur_lay)[1:nb_fact_lay] df_cur_lay <- cur_lay$Df[nrow(cur_lay)] SS_cur_lay <- sum(cur_lay$"Sum Sq") for (ifact in 1:nb_fact_lay) { df_cur_fact <- cur_lay$Df[ifact] SS_cur_fact <- cur_lay$"Sum Sq"[ifact] Fv_cur_fact <- round(cur_lay$"F value"[ifact],digits=2) pv_cur_fact <- cur_lay$"Pr(>F)"[ifact] pvseuils <- c(Inf, 0.1,0.05,0.01,0.001) pv_thr <- last(pvseuils[pv_cur_fact < pvseuils]) etasquare <- round(SS_cur_fact/SS_cur_lay,digits=3) petasquare <- round(SS_cur_fact/(SS_cur_fact+SS_cur_lay),digits=3) if (maineff!="off") { print(paste(nm_fact_lay[ifact]," (F(", df_cur_fact, ",", df_cur_lay, ")=", Fv_cur_fact, ", p=", signif(pv_cur_fact,3), "<", pv_thr, ", eta-square=", etasquare, ", eta-square=", petasquare, ")", sep="")) } } } if (maineff=="off") { print(paste(nm_fact_lay[ifact]," (F(", df_cur_fact, ",", df_cur_lay, ")=", Fv_cur_fact, ", p=", signif(pv_cur_fact,3), "<", pv_thr, ", eta-square=", etasquare, ", peta-square=", petasquare, ")", sep="")) } #This line was in the loops (two lines above here, right after "etasquare <- ..."). This would report main effects. } # BOOTSTRAP require(MASS) ttestboot <- function (a, b, N=1000, paired=T, titre="") { a <- na.omit(a) b <- na.omit(b) if (paired==T) { dist <- c(a-b, b-a) samples <- c() for (i in 1:N) { s <- sample(dist, size=length(a), replace=TRUE) samples[i] <- t.test(s, mu=0)$statistic } t0 <- t.test(a-b, mu=0)$statistic } else { dist <- c(a, b) samples <- c() for (i in 1:N) { s1 <- sample(dist, size=length(a), replace=TRUE) s2 <- sample(dist, size=length(a), replace=TRUE) samples[i] <- t.test(s1, s2)$statistic } t0 <- t.test(a, b)$statistic } p1 <- min(sum(t0samples))/length(samples) p2 <- sum(abs(t0)100] <- 100 dist <- data.frame ( rbind(dist, dist2) ) # for (x in 1:length(table(factor(g)))) { # gL <- names(table(factor(g)))[x] # # dist2 <- data.frame(a=a, f=f, g=g, s=paste(s, gL, sep="-") ) # dist2$f[dist2$g==gL & dist2$f==fL1] <- NA # dist2$f[dist2$g==gL & dist2$f==fL2] <- fL1 # dist2$f[is.na(dist2$f)] <- fL2 # dist <- data.frame ( rbind(dist, dist2) ) # } samples <- c() for (i in 1:N) { ss <- sample(names(table(factor(dist$s))), size=Ns, replace=TRUE) T <- subset(dist, s %in% ss) an <- aov(a ~ factor(f) * factor(g) + Error(factor(s)/(factor(f)*factor(g))) , data=T) lay <- eval(parse(text=paste("summary(an)$ \"",names(summary(an)[4]),"\"",sep="")))[[1]] samples[i] <- lay$"F value"[1] } T0 <- subset(dist, s %in% names(table(factor(dist0$s)))) an <- aov(a ~ f * g + Error(s/(f*g)) , data=T0) lay <- eval(parse(text=paste("summary(an)$ \"",names(summary(an)[4]),"\"",sep="")))[[1]] F0 <- lay$"F value"[1] p <- sum(F0