#Analyses using MINT
library(evolqg)
library(geomorph)

#simulate data with known covariance
sim.data<-function(n,s){
  a<-matrix(rnorm(n*nrow(s)),ncol=nrow(s))
  b<-a%*%chol(s)
  return(b)
}

#builds with Sw and Sb with random noise in all components
s.build<-function(p,sb,sw,mod.gp){
  s<-matrix(sb,nrow=p, ncol=p)
  for (i in 1:nlevels(mod.gp)){
    s[which(mod.gp==levels(mod.gp)[i]),which(mod.gp==levels(mod.gp)[i])]<-sw
  }
  low<-s[lower.tri(s)]
  for (i in 1:length(low)){
    low[i]<-rnorm(1,mean=low[i],sd=0.01)  
  }
  s[lower.tri(s)]<-low
  s[upper.tri(s)]<-t(s)[upper.tri(s)] 
  diag(s)<-rnorm(p,mean=1,sd=0.01)
  s<-matrix(Matrix::nearPD(s)$mat,nrow=p,ncol=p)
  return(s)
}

#FOR random covariance matrix: 2 module
N <- c(100)
p <- c(16,32)
nsimul <- 100
###### Power Simulation
s.b <- 0.6
s.w <- c(0,.025,.05,.075,.1) #NOTE, this is ADDED to SB

MINT.Power <- MINT.Misspec <- array(NA, dim = c(length(p), length(N),length(s.w)))
row.names(MINT.Power) <- row.names(MINT.Misspec) <- paste("p_",p,sep="")
colnames(MINT.Power) <- colnames(MINT.Misspec) <- paste("N_",N,sep="")
dimnames(MINT.Power)[[3]] <- dimnames(MINT.Misspec)[[3]] <- paste("SW_",s.w,sep="")
for (i in 1:length(p)){
  for (j in 1:length(N)){
    #hypotheses
    m2 <- gl(2,p[i]/2)
    m4 <- gl(4,p[i]/4)
    m8 <- gl(8,p[i]/8)
    H.modules <- cbind(model.matrix(~m2-1),model.matrix(~m4-1),model.matrix(~m8-1))
    
    #data 
    for (ii in 1:length(s.w)){
      S.in <- s.build(p[i],sb=s.b,sw=(s.b+s.w[ii]),mod.gp=m4)
      Y <- lapply(1:nsimul, function(k) sim.data(n=N[j],S.in))
      MINT.out <- lapply(1:nsimul, function(k)  JackKnifeMINT(Y[[k]],H.modules, n = N) )
      
      #best model: Type 1 is false positives , 'power' = correct model, misspec = incorrect model  
      best <- unlist(lapply(1:nsimul, function(k) MINT.out[[k]]$.id[1])) #name of best model
      
      MINT.Power[i,j,ii] <- ifelse(ii==1,  1-length(which(best== "null")) / nsimul ,
                                   length(which(best=="m41_m42_m43_m44")) /nsimul ) 
      MINT.Misspec[i,j,ii] <- 1-length(which(best=="m41_m42_m43_m44")) /nsimul
      
     
      print(paste("i =",i))
      print(paste("j =",j))
      print(paste("ii =",ii))
    }
  }
}

write.csv(MINT.Power,"mint.power-4gp-.6sb.csv")
write.csv(MINT.Misspec,"mint.misspec-4gp-.6sb.csv")


