## Demonstration
require(OpenMx)
source("GW-SEM.R")


# Run the Ordinal 2 factor GW-SEM
start <- Sys.time()
snpCovs(FacModelData = "Ord2fac.txt", vars = c("y1","y2","y3","y4","x1","x2","x3","x4"), controls = c("c1","c2","c3"), SNPdata = "SNP.txt", output = "Ord2fac")
twofacDWLS(itemData = "Ord2fac.txt", snpCov = "Ord2faccov.txt", snpWei = "Ord2facwei.txt", f1Names = c("y1","y2","y3","y4"), f2Names = c("x1","x2","x3","x4"), Controls = c("c1","c2","c3"), output = "Ord2FacOut.txt")	
end <- Sys.time()
end-start



twoFac <- read.table("Ord2FacOut.txt")
par(mfrow = c(2,1), mar = c(0, 5, 0, 0) + 0.1)
	plot(-log10(twoFac$snp2fac1.p), ylim = c(0, 17), ylab = "Factor 1  \n -Log10 (P value)", 
	     main = " ", xlab = " ", xaxt = "n", pch = 16)
	abline(h = 5, col = "blue", lwd = 2)
	abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)

	plot(-log10(twoFac$snp2fac2.p), ylim = c(0, 17), ylab = "Factor 2  \n -Log10 (P value)", 
	     main = " ", xlab = " ", xaxt = "n", pch = 16)
	abline(h = 5, col = "blue", lwd = 2)
	abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)
	
	
	
head(twoFac[order(twoFac$snp2fac1.p),c("snp2fac1","snp2fac2","snp2fac1SE","snp2fac2SE","snp2fac1.p","snp2fac2.p")], 10)
head(twoFac[order(twoFac$snp2fac2.p),c("snp2fac1","snp2fac2","snp2fac1SE","snp2fac2SE","snp2fac1.p","snp2fac2.p")], 10)
	
	
# Examine the Estimates from the one factor GW-SEM and the residuals GW-SEM
twoFacEst  <- read.table("twoFacEstimates.txt")
twoFacEst[c("rs2178","rs2014","rs6001"),]



           # Make a QQ plot for the 2 Factor GW-SEM
           require(car)
           twoFacRes <- read.table("Ord2FacOut.txt")
           par(mfrow= c(1,2))
           qqPlot(twoFacRes$snp2fac1.t, grid = F, main = "QQ-Plot for Factor 1", 
                  ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")
           qqPlot(twoFacRes$snp2fac2.t, grid = F, main = "QQ-Plot for Factor 2", 
                  ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")