# Type I error and power simulations: old vs. new RR methods univariate & multivariate library(ape) library(phytools) library(nlme) library(geiger) library(geomorph) nsim = 1000 n=50 a <- 1; a.pos <- 1.15 phy <- pbtree(n = n) x <- lapply(1:nsim, function(j) fastBM(tree = phy,mu = 4,a = 0)) y <- lapply(1:nsim, function(j) a*x[[j]] + rnorm(n, 0, .1)) #TYPE 1 error test y.pos <- lapply(1:nsim, function(j) (a.pos)*x[[j]] + rnorm(n, 0, .1)) #RR: power df <- lapply(1:nsim, function(j) data.frame(x[[j]], y[[j]], y.pos[[j]])) ## Standard method fit <- lapply(1:nsim, function(j) gls(y..j..~x..j.., correlation=corBrownian(phy=phy), data=df[[j]])) T <- lapply(1:nsim, function(j) (coef(fit[[j]])[2]-1) / summary(fit[[j]])$tTable[2,2]) prob.T <- unlist(lapply(1:nsim, function(j) dt(T[[j]],df = (nrow(df[[j]])-1)))) length(which(prob.T< 0.05))/nsim # Type 1 ~ 0.047 fit.RR <- lapply(1:nsim, function(j) gls(y.pos..j..~x..j.., correlation=corBrownian(phy=phy), data=df[[j]])) T.RR <- lapply(1:nsim, function(j) (coef(fit.RR[[j]])[2]-1) / summary(fit.RR[[j]])$tTable[2,2]) prob.T.RR <- unlist(lapply(1:nsim, function(j) dt(T.RR[[j]],df = (nrow(df[[j]])-1)))) length(which(prob.T.RR< 0.05))/nsim # Power ~ 0.931 #SDDI Method mn.sz <- lapply(1:nsim, function(j) rowMeans(cbind(y[[j]],x[[j]]))) rel.x <- lapply(1:nsim, function(j) x[[j]] / mn.sz[[1]] ) rel.y <- lapply(1:nsim, function(j) y[[j]] / mn.sz[[1]] ) SDDI <- lapply(1:nsim, function(j) apply(cbind((rel.y[[j]]),(rel.x[[j]])),1,dist)) for (j in 1:nsim){ SDDI[[j]][which(x[[j]] > y[[j]])] <- -1*SDDI[[j]][which(x[[j]] > y[[j]])] } fit.SDDI <- lapply(1:nsim, function(j) gls(SDDI~mn.sz, correlation = corBrownian(phy=phy), data = data.frame(SDDI=SDDI[[j]], mn.sz = mn.sz[[j]]))) P.SDDI <- unlist(lapply(1:nsim, function(j) anova(fit.SDDI[[j]])[2,3])) length(which(P.SDDI< 0.05))/nsim # Type 1 ~ 0.057 #RR power mn.sz.RR <- lapply(1:nsim, function(j) rowMeans(cbind(y.pos[[j]],x[[j]]))) rel.x.RR <- lapply(1:nsim, function(j) x[[j]] / mn.sz.RR[[1]] ) rel.y.RR <- lapply(1:nsim, function(j) y.pos[[j]] / mn.sz.RR[[1]] ) SDDI.RR <- lapply(1:nsim, function(j) apply(cbind((rel.y.RR[[j]]),(rel.x.RR[[j]])),1,dist)) for (j in 1:nsim){ SDDI.RR[[j]][which(x[[j]] > y.pos[[j]])] <- -1*SDDI.RR[[j]][which(x[[j]] > y.pos[[j]])] } fit.SDDI.RR <- lapply(1:nsim, function(j) gls(SDDI.RR~mn.sz.RR, correlation = corBrownian(phy=phy), data = data.frame(SDDI.RR=SDDI.RR[[j]], mn.sz.RR = mn.sz.RR[[j]]))) P.SDDI.RR <- unlist(lapply(1:nsim, function(j) anova(fit.SDDI.RR[[j]])[2,3])) length(which(P.SDDI.RR< 0.05))/nsim # Power ~ 0.97 ################ Multivariate: p=5 #DONE AS MOSIMANN SHAPE VARIABLES nsim = 1000 n=50 phy <- pbtree(n = n) p <- 5 R <- matrix(0.7,nrow=p+1,ncol=p+1); diag(R) <- 1 x.o <- lapply(1:nsim, function(j) sim.char(phy = phy, par = R,root=10)[,,1]) mn.sz <- lapply(1:nsim, function(j) x.o[[j]][,1]) x <- lapply(1:nsim, function(j) x.o[[j]][,-1]/2) fit <- lapply(1:nsim, function(j) lm(x[[j]]~mn.sz[[j]]) ) y <- lapply(1:nsim, function(j) x[[j]] + mvtnorm::rmvnorm(n,sigma= R[2:(p+1),2:(p+1)]/10) ) a <-lapply(1:nsim, function(j) 1.15*fit[[j]]$coefficients) y.pos <- lapply(1:nsim, function(j) model.matrix(~mn.sz[[j]])%*%a[[j]] + fit[[j]]$residuals +mvtnorm::rmvnorm(n,sigma= R[2:(p+1),2:(p+1)]/100000)) #NOTE: random error added in scaled proportion to input signal (so not to swamp) x <- lapply(1:nsim, function(j) x[[j]] / mn.sz[[j]]) y <- lapply(1:nsim, function(j) y[[j]] / mn.sz[[j]]) y.pos <- lapply(1:nsim, function(j) y.pos[[j]] / mn.sz[[j]]) for (j in 1:nsim){rownames(y[[j]]) <- rownames(y.pos[[j]]) <- rownames(x[[j]])} df.s <- lapply(1:nsim, function(j) list(x[[j]], y[[j]], y.pos[[j]])) #SDDI Method multi ShDDI <- ShD.dist <- lapply(1:nsim, function(j) sqrt(apply((df.s[[j]][[1]]-df.s[[j]][[2]])^2,1,sum)) ) #PLS method PLS <- lapply(1:nsim, function(j) two.b.pls(A1 = x[[j]], A2 = y[[j]], print.progress = FALSE, iter=1)) res <- lapply(1:nsim, function(j) cbind(PLS[[j]]$XScores[,1],PLS[[j]]$YScores[,1])) flip <- lapply(1:nsim, function(j) which(apply(res[[j]],1,which.min)==1)) for (j in 1:nsim){ ShD.dist[[j]][flip[[j]]] <- -1*ShD.dist[[j]][flip[[j]]] } fit.ShD <- lapply(1:nsim, function(j) gls(ShD~mn.sz, correlation = corBrownian(phy=phy), data = data.frame(ShD=ShD.dist[[j]], mn.sz = mn.sz[[j]]))) P.ShD <- unlist(lapply(1:nsim, function(j) anova(fit.ShD[[j]])[2,3])) length(which(P.ShD< 0.05))/nsim # Type 1 ~ 0.054 ### NOTE: if one does not flip, type I error remains the same #RR: SDDI Method multi power ShD.dist.RR <- lapply(1:nsim, function(j) sqrt(apply((df.s[[j]][[1]]-df.s[[j]][[3]])^2,1,sum)) ) PLS.RR <- lapply(1:nsim, function(j) two.b.pls(A1 = x[[j]], A2 = y.pos[[j]], print.progress = FALSE, iter=1)) res.RR <- lapply(1:nsim, function(j) cbind(PLS.RR[[j]]$XScores[,1],PLS.RR[[j]]$YScores[,1])) flip.RR <- lapply(1:nsim, function(j) which(apply(res.RR[[j]],1,which.min)==1)) for (j in 1:nsim){ ShD.dist.RR[[j]][flip.RR[[j]]] <- -1*ShD.dist.RR[[j]][flip.RR[[j]]] } fit.ShD.RR <- lapply(1:nsim, function(j) gls(ShD~mn.sz, correlation = corBrownian(phy=phy), data = data.frame(ShD=ShD.dist.RR[[j]], mn.sz = mn.sz[[j]]))) P.ShD.RR <- unlist(lapply(1:nsim, function(j) anova(fit.ShD.RR[[j]])[2,3])) length(which(P.ShD.RR< 0.05))/nsim # Power ~ 0.98 #NOTE: if one does not flip, power remains the same ########### #PLOT Type1.Std <- length(which(prob.T< 0.05))/nsim Type1.SDDI <- length(which(P.SDDI< 0.05))/nsim Type1.multi <- length(which(P.ShD< 0.05))/nsim Power.Std <- length(which(prob.T.RR< 0.05))/nsim Power.SDDI <- length(which(P.SDDI.RR< 0.05))/nsim Power.multi <- length(which(P.ShD.RR< 0.05))/nsim sim.val <- cbind(Type1.Std,Type1.SDDI,Type1.multi,Power.Std,Power.SDDI,Power.multi) barplot(sim.val, ylim = c(0,1)) abline(h=0.05)