--- title: "R Notebook" output: html_document: df_print: paged --- ```{r} Xs = test2 nperm=100 ``` ```{r} ix = c(4,3,6,5,8,7) iy = 1 r = corls(Xs,permutations = how(nperm = nperm)) rx = r[ix,ix] ry = r[iy,ix] rxi = r[ix[(ix %% 2)>0],ix[(ix %% 2)>0]] ryi = r[iy,ix[(ix %% 2)>0]] ``` ```{r} r2 = (solve(rx) %*% ry) * ry r2 ``` ```{r} sum(r2) ``` ```{r} r2 = (solve(rxi) %*% ryi) * ryi r2 ``` ```{r} sum(r2) ``` ```{r} v = varls(Xs,permutations = how(nperm = nperm)) vt = v[c(iy,ix[(ix %% 2)>0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1) vr = v[c(iy+1,ix[(ix %% 2)==0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1) vt vr ``` ```{r} #diag(vr)=0 vp = vt - vr vp ``` ```{r} d = sqrt(diag(vp)) n = d %o% d r = vp / n r ``` ```{r} rx = r[ix[(ix %% 2)==0]/2,ix[(ix %% 2)==0]/2] ry = r[iy,ix[(ix %% 2)==0]/2] r2 = (solve(rx) %*% ry) * ry r2 ``` ```{r} sum(r2) ``` ```{r} rp = solve(r) rp = rp / sqrt(diag(rp) %o% diag(rp)) rp rp^2 ``` ```{r} library(ProcMod) r2s = seq(from=0.0, to = 0.9, by = 0.1) n.sim=1000 r2.obs = matrix(0,nrow = n.sim,ncol = length(r2s)) d1 = simulate_matrix(10,2) for (j in seq_along(r2s)) { #print(c("r2",as.character(r2s[j]))) for (i in seq_len(n.sim)) { # print(c("round",as.character(i))) d2 = simulate_correlation(d1,1,r2s[j]) #r2.obs[i,j] = cor5(d1,d2)^2 r2.obs[i,j] = ((max(0,.covls(d1,d2)-.rcovls(d1,d2)))/(sqrt(.covls(d1,d1)-.rcovls(d1,d1))*sqrt(.covls(d2,d2)-.rcovls(d2,d2))))^2 # v = varls(d1,d2,permutations = 10) # s = sqrt(diag(v)) # r2.obs[i,j] = (v / outer(s,s))[1,2]^2 } } boxplot(r2.obs) cbind(R2_theo=r2s, R2_obs=colMeans(r2.obs), R2_sd=apply(r2.obs,MARGIN = 2,FUN = sd)) ``` ```{r} a= simulate_matrix(40,1) b= simulate_matrix(40,2) x = varls(a,b,permutations = NULL) y = corls(a,b,permutations = 100) y^2 m = x - y #(sqrt(outer(diag(x),diag(x)))/m)^2 #sqrtm(outer(diag(x),diag(x))) ``` 5 -> 6 10 -> 12 6 20 -> 28 14 30 -> 44 22 40 -> 60 30 20 ```{r} i=c(5,10,20,30,40) j=c(6,12,28,44,60) lm(j~i) plot(i,j,xlim=c(0,50),ylim = c(-10,70)) abline(lm(j~i)) plot(sqrt(i),sqrt(j)) lm(sqrt(j)~sqrt(i)) plot(i^2,j^2) lm(i^2~j^2) ``` ```{r} a = simulate_matrix(20,1) b = vector(mode = "list",9) for (i in seq_along(b)) { b[[i]] = simulate_correlation(a,1,i/10) } b = c(list(a),b) var(do.call(cbind,b)) var2(do.call(cbind,b)) v.perm=varls(as.procmod.frame(b),permutations = 100,rcovls = TRUE) v.perm m = apply(attributes(v.perm)\$rcov,MARGIN = c(1,2),mean,na.rm=TRUE)^2 mean(m[upper.tri(m,diag = TRUE)],na.rm=TRUE) ``` ```{r} library(ProcMod) ps = c(seq(1,10,by=1),seq(15,100,by=5),1000) qs = ps ns = seq(5,20, by=5) n.sim=1000 r2.obs = array(0,dim = c(length(ns),length(ps),length(qs),n.sim)) for (ni in seq_along(ns)) { N=ns[ni] for (pi in seq_along(ps)) { P = ps[pi] for (qi in seq_along(qs)) { Q = qs[qi] for (si in 1:n.sim) { not_ok = TRUE while(not_ok) { d1 = simulate_matrix(N,P,equal.var = TRUE) d2 = simulate_matrix(N,Q,equal.var = TRUE) r2=tryCatch({v = varls(d1,d2,permutations = 0) s = sqrt(diag(v)) (v / outer(s,s))[1,2]^2 }, error=function(cond) {#message("Redo simulation") NA} ) r2.obs[ni,pi,qi,si] = r2 not_ok=is.na(r2) } } } } } ``` ```{r} R2_obs=apply(r2.obs,MARGIN = c(1,2,3),FUN = mean) R2_sd=apply(r2.obs,MARGIN = c(1,2,3),FUN = sd) dimnames(R2_obs) = list(ns,ps,qs) N = rep(ns,length(ps)*length(qs)) P = rep(rep(ps,rep(length(ns),length(ps))),length(qs)) Q = rep(qs,rep(length(ns)*length(ps),length(qs))) R2_obs.tab = tibble(N=N, P=P, Q=Q, Pmin = mapply(min,P,Q), Qmax = mapply(max,P,Q), mean=as.numeric(R2_obs), sd=as.numeric(R2_sd)) R2_obs.tab\$NPmin = mapply(min,R2_obs.tab\$Pmin,N) R2_obs.tab\$NPmax = mapply(max,R2_obs.tab\$Pmin,N) write_csv2(R2_obs.tab,"R2_H0.csv") #p = rep(ps,rep(length(qs),length(ps))) #q = rep(qs,length(ps)) #R2_obs.tn = data.frame(p,q,N=N,mean=as.numeric(R2_obs),sd=as.numeric(R2_sd)) ``` ```{r} library(expm) ``` # Modèle 1 : l'axe suivant se prend le reste ```{r} R2.alea.predict.single = function(N,nP,nQ) { # N : nombre de points # nP : Nb de var dans petit tableau # nQ : Nb de var dans grand tableau # Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1) Axes.main = min(nP,(N-1)) i = 0:(Axes.main-1) Q = 1/(N-1) P = (N-2)/(N-1) Pis = P**i R2.main = sum(Pis)*Q # Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ; Axes.comp = max(nP,nQ) - Axes.main if (Axes.comp >0) { i = 0:(Axes.comp-1) Q = 1/(N-1)^2 P = 1 - Q Pis = (P**i) R2.main = P**Axes.comp*R2.main+Q*sum(Pis) } return(R2.main) } ``` # Modèle 2 : l'axe suivant se prend le reste ```{r} R2.alea.predict.single = function(N,nP,nQ) { # N : nombre de points # nP : Nb de var dans petit tableau # nQ : Nb de var dans grand tableau # Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1) Axes.main = min(nP,(N-1)) Rres = numeric(Axes.main) R2.main=0 for (i in 1:Axes.main) { Rres[i]=1-R2.main if (i==1) k = 0 else k = sum(Rres[1:(i-1)]^2) # print(k / (N-1)^2 * Rres[i]) R2.main = R2.main + Rres[i] * ( 1/(N-1) - k/(N-1)^2) } # R2.main = R2.main - 2*Axes.main / (N-1)^2 # print(2*Axes.main / (N-1)^2) # Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ; k=sum(Rres^2) + (1 - R2.main)^2 Axes.comp = max(nP,nQ) - Axes.main if (Axes.comp >0) { for (i in 1:Axes.comp) { #print(R2.main) R2.main = R2.main + (1-R2.main) * k/(N-1)^2 } } return(R2.main) } R2.alea.predict =Vectorize(R2.alea.predict.single, vectorize.args = c("N","nP","nQ")) ``` # Modèle 3 : l'axe suivant se prend le reste ```{r} R2.alea.predict.single = function(N,nP,nQ) { # N : nombre de points # nP : Nb de var dans petit tableau # nQ : Nb de var dans grand tableau # Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1) Axes.main = min(nP,(N-1)) Rres = numeric(Axes.main) R2.main=0 for (i in 1:Axes.main) { Rres[i]=1-R2.main R2.main = R2.main + Rres[i] * ( 1/(N-1) - (Axes.main - i )/(N-1)^2) } # R2.main = R2.main - 2*Axes.main / (N-1)^2 # print(2*Axes.main / (N-1)^2) # Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ; k=sum(Rres) + (1 - R2.main) Axes.comp = max(nP,nQ) - Axes.main if (Axes.comp >0) { for (i in 1:Axes.comp) { #print(R2.main) R2.main = R2.main + (1-R2.main) * k/(N-1)^2 } } return(R2.main) } R2.alea.predict =Vectorize(R2.alea.predict.single, vectorize.args = c("N","nP","nQ")) ``` ```{r} library(tidyverse) R2_obs.tab = read_csv2("R2_H0.csv") ``` ```{r} R2_obs.tab.PQ = R2_obs.tab %>% filter(Qmax==Pmin) ``` ```{r} with(data = R2_obs.tab.PQ, plot(mean, R2.alea.predict(N,Pmin,Qmax), col=rainbow(4)[N/5], pch=(1+(Qmax==Pmin)), cex=0.5, log="") ) abline(b=1,a=0) abline(h=1,lty=2) abline(v=1,lty=2) abline(lm(R2.alea.predict(N,Pmin,Qmax)~mean,data = R2_obs.tab.PQ)) #summary(lm(mean~R2.alea.predict(N,Pmin,Qmax),data = R2_obs.tab.PQ)) legend("bottomright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5) ``` ```{r} with(data = R2_obs.tab.PQ, plot(mean, R2.alea.predict(N,Pmin,Qmax)-mean, col=rainbow(4)[N/5], pch=(1+(Qmax==Pmin)), cex=0.5, log="") ) abline(h=1) abline(v=2) abline(v=c(5,10,15,20)*2.5,lty=2) #abline(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ)) legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5) #summary(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ)) ``` ```{r} with(data = R2_obs.tab.PQ, plot(Pmin, R2.alea.predict(N,Pmin,Qmax)-mean, col=rainbow(4)[N/5], pch=(1+(Qmax==Pmin)), cex=0.5, log="x") ) abline(h=1) abline(v=2) abline(v=c(5,10,15,20)*4,lty=2,col=rainbow(4)) abline(v=c(5,10,15,20),lty=3,col=rainbow(4)) #abline(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ)) legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5) #summary(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ)) ``` ```{r} with(data = R2_obs.tab.PQ, plot(Pmin, mean, col=rainbow(4)[N/5], pch=(1+(Qmax==Pmin)), cex=0.5, log="x") ) abline(h=1) legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5) ```