A bit more

.Rbuildignore

@@ -1,2 +1,3 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^\.travis\.yml$

.gitignore

@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 .Ruserdata
+/NAMESPACE

.travis.yml

@@ -0,0 +1,40 @@
+# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r
+
+language: R
+cache: packages
+dist: xenial
+latex: true
+
+addons:
+  apt:
+    sources:
+    - ubuntu-toolchain-r-test
+    packages:
+    - g++-8
+
+env:
+  global:
+    - MAKEFLAGS="-j 2"
+
+before_install:
+  - mkdir -p ~/.R/
+  - echo "CXX14 = g++-8 -fPIC -flto=2" >> ~/.R/Makevars
+  - echo "CXX14FLAGS = -mtune=native -march=native -Wno-unused-variable -Wno-unused-function -Wno-unused-local-typedefs -Wno-ignored-attributes -Wno-deprecated-declarations -Wno-attributes -O3" >> ~/.R/Makevars
+
+install:
+- ./tools/travis/install-$TRAVIS_OS_NAME
+
+before_script:
+- export _R_CHECK_FORCE_SUGGESTS_=false
+- echo "pkgbuild::compile_dll(); devtools::document(roclets = c('rd', 'collate', 'namespace'))" | R --vanilla
+
+script: 
+  - |
+    R CMD build --no-build-vignettes .
+    travis_wait 59 R CMD check --ignore-vignettes gwsem*tar.gz
+
+after_success:
+  - travis_wait 40 Rscript -e 'covr::codecov()'
+
+after_failure:
+  - cat gwsem.Rcheck/00*

DESCRIPTION

@@ -10,4 +10,10 @@ License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true
 Language: en-US
-Imports: OpenMx
+LinkingTo: BH (>= 1.69.0-1)
+Depends: 
+    OpenMx (>= 2.13.2.187)
+Suggests: 
+    testthat
+RoxygenNote: 6.1.1
+SystemRequirements: GNU make

R/model.R

@@ -0,0 +1,245 @@
+require(OpenMx)
+require(MASS)
+require(psych)
+require(data.table)
+
+
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+
+ oneFacGWAS <- function(phenoData, snpData, itemNames, covariates = NULL, nSNP = NA, fitfun = "WLS", minMAF = .01, snpFileType = 'bgen', out = "out"){
+ 	minVar <- 2*minMAF*(1-minMAF)
+
+    fac <- matrix(1, length(itemNames))
+    thr <- matrix(1, length(itemNames))
+
+    for(i in 1:length(itemNames)){
+   fac[i] <-  is.factor(phenoData[,itemNames][,i])
+   thr[i] <-  nlevels(phenoData[,itemNames][,i])-1
+    }
+    thr[thr< 0] <- 0
+    maxThr <- max(thr)
+
+phenoData$snp <- rbinom(dim(phenoData)[1], 2, .5) # create placeholder
+latents   <- c("F")
+lambda    <- mxPath(from=latents, to=itemNames,values=1, free = T, labels = paste("lambda", itemNames, sep = "_")  )
+snpMu     <- mxPath(from = "one", to = "snp" , labels = "snpMean")
+snpBeta   <- mxPath(from = "snp", to = "F", labels = "snpReg", values = 0, free = T)
+snpres    <- mxPath(from = "snp", arrows=2, values=1, free = T, labels = paste("snp", "res", sep = "_"))
+resid     <- mxPath(from = c(itemNames), arrows=2, values=1, free = c(fac==0), labels = paste(c(itemNames), "res", sep = "_"))
+facRes    <- mxPath(from=latents, arrows=2,free=F, values=1.0, labels = "facRes")
+covMean   <- mxPath(from = "one", to = covariates, free=FALSE, labels = paste0('data.',covariates)) 
+cov2item  <- mxPath(from = covariates, to = itemNames, connect = "all.pairs", labels = paste(rep(covariates, each = length(itemNames)), itemNames, sep = "_2_"))
+itemMean  <- mxPath(from = 'one', to = itemNames, free= c(fac==0), values = 0, labels = paste0(itemNames, "Mean"))
+
+if(maxThr>0) thresh    <- mxThreshold(itemNames[c(fac==1)], nThresh=c(thr[fac==1]), free = T , labels = paste(rep(itemNames[c(fac==1)], each = maxThr), "_Thr_", 1:maxThr, sep = ""), values=mxNormalQuantiles(1))
+
+
+dat       <- mxData(observed=phenoData, type="raw")                                              ## options for min variance/MAF
+
+if(maxThr==0 & fitfun == "WLS") fun      <- mxFitFunctionWLS(allContinuousMethod= "marginals")
+if(maxThr>0 & fitfun == "WLS") fun       <- mxFitFunctionWLS()
+if(fitfun == "FIML") fun                 <- mxFitFunctionML()
+
+
+ oneFacPre <- mxModel("OneFac", type="RAM",
+       manifestVars = c("snp", itemNames),
+       latentVars = c(latents, covariates),
+ 	  lambda, snpMu, snpBeta, snpres, resid, facRes, covMean, cov2item,
+ 	  itemMean, dat, fun  )
+
+
+if(maxThr>0) oneFacPre <- mxModel(oneFacPre, name = "OneFac", thresh  )
+
+ if(snpFileType == "bgen"){
+ Compute <- mxComputeLoop(list(
+ 	 mxComputeSetOriginalStarts(),
+      mxComputeLoadData("OneFac", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+ 	 mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+                                ), i=1:nSNP)}
+
+ if(snpFileType == "pgen"){
+ Compute <- mxComputeLoop(list(
+      mxComputeSetOriginalStarts(),
+      mxComputeLoadData("OneFac", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+      mxComputeLoadContext(path=paste(snpData, "pvar", sep = "."),column=1:3, sep='\t'),
+      mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+                                ), i=1:nSNP)
+  }
+
+
+ oneFac <- mxModel(oneFacPre, name = "OneFac", Compute  )
+
+ oneFacFit <- mxRun(oneFac)
+ summary(oneFacFit)
+
+ }
+
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+
+
+
+
+
+ oneFacResGWAS <- function(phenoData, snpData, itemNames , factor = F, res = itemNames, covariates = NULL, nSNP = NA, fitfun = "WLS", minMAF = .01, snpFileType = 'bgen', out = "res"){
+ 	minVar <- 2*minMAF*(1-minMAF)
+
+    fac <- matrix(1, length(itemNames))
+    thr <- matrix(1, length(itemNames))
+
+    for(i in 1:length(itemNames)){
+   fac[i] <-  is.factor(phenoData[,itemNames][,i])
+   thr[i] <-  nlevels(phenoData[,itemNames][,i])-1
+    }
+    thr[thr< 0] <- 0
+    maxThr <- max(thr)
+
+
+phenoData$snp <- rbinom(dim(phenoData)[1], 2, .5) # create placeholder
+latents   <- c("F")
+lambda    <- mxPath(from=latents, to=itemNames,values=1, free = T, labels = paste("lambda", itemNames, sep = "_")  )
+snpMu     <- mxPath(from = "one", to = "snp" , labels = "snpMean")
+snpFac  <- mxPath(from = "snp", to = "F", labels = "snp2fac", values = 0, free = factor)
+snpItemRes   <- mxPath(from = "snp", to = res, labels = paste("snp", res, sep = "2"), values = 0, free = T)
+snpres    <- mxPath(from = "snp", arrows=2, values=1, free = T, labels = paste("snp", "res", sep = "_"))
+resid     <- mxPath(from = c(itemNames), arrows=2, values=1, free = c(fac==0), labels = paste(c(itemNames), "res", sep = "_"))
+facRes    <- mxPath(from=latents, arrows=2,free=F, values=1.0, labels = "facRes")
+covMean   <- mxPath(from = "one", to = covariates, free=FALSE, labels = paste0('data.',covariates)) 
+cov2item  <- mxPath(from = covariates, to = itemNames, connect = "all.pairs", labels = paste(rep(covariates, each = length(itemNames)), itemNames, sep = "_2_"))
+itemMean  <- mxPath(from = 'one', to = itemNames, free= c(fac==0), values = 0, labels = paste0(itemNames, "Mean"))
+
+
+if(maxThr>0) thresh    <- mxThreshold(itemNames[c(fac==1)], nThresh=c(thr[fac==1]), free = T , labels = paste(rep(itemNames[c(fac==1)], each = maxThr), "_Thr_", 1:maxThr, sep = ""), values=mxNormalQuantiles(1))
+dat       <- mxData(observed=phenoData, type="raw")                                              ## options for min variance/MAF
+
+if(maxThr==0 & fitfun == "WLS") fun      <- mxFitFunctionWLS(allContinuousMethod= "marginals")
+if(maxThr>0 & fitfun == "WLS") fun       <- mxFitFunctionWLS()
+if(fitfun == "FIML") fun                 <- mxFitFunctionML()
+
+
+ oneFacPre <- mxModel("OneFacRes", type="RAM",
+   manifestVars = c("snp", itemNames),
+   latentVars = c(latents, covariates),
+   lambda, snpMu, snpFac, snpItemRes, snpres, resid, facRes, covMean, cov2item,
+   itemMean, dat, fun  )
+
+if(maxThr>0) oneFacPre <- mxModel(oneFacPre, name = "OneFacRes", thresh  )
+
+	 if(snpFileType == "bgen"){
+	 Compute <- mxComputeLoop(list(
+	 	 mxComputeSetOriginalStarts(),
+	      mxComputeLoadData("OneFacRes", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+	 	 mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+	      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+	                                ), i=1:nSNP)}
+
+	 if(snpFileType == "pgen"){
+	 Compute <- mxComputeLoop(list(
+	      mxComputeSetOriginalStarts(),
+	      mxComputeLoadData("OneFacRes", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+	      mxComputeLoadContext(path=paste(snpData, "pvar", sep = "."),column=1:3, sep='\t'),
+	      mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+	      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+	                                ), i=1:nSNP)  }
+
+ oneFac <- mxModel(oneFacPre, name = "OneFacRes", Compute  )
+ oneFacFit <- mxRun(oneFac)
+ summary(oneFacFit)
+
+	 }
+
+
+
+
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+
+
+twoFacGWAS <- function(phenoData, snpData, F1itemNames, F2itemNames, covariates = NULL, nSNP = NA, fitfun = "WLS", minMAF = .01, snpFileType = 'bgen', out = "out"){
+ 	minVar <- 2*minMAF*(1-minMAF)
+
+    fac <- matrix(1, length(c(F1itemNames,F2itemNames)))
+    thr <- matrix(1, length(c(F1itemNames,F2itemNames)))
+
+    for(i in 1:length(c(F1itemNames,F2itemNames))){
+   fac[i] <-  is.factor(phenoData[,c(F1itemNames,F2itemNames)][,i])
+   thr[i] <-  nlevels(phenoData[,c(F1itemNames,F2itemNames)][,i])-1
+    }
+    thr[thr< 0] <- 0
+    maxThr <- max(thr)
+
+phenoData$snp <- rbinom(dim(phenoData)[1], 2, .5) # create placeholder
+latents   <- c("F1", "F2")
+lambda1    <- mxPath(from="F1", to=F1itemNames,values=1, labels = paste("lambda", F1itemNames, sep = "_")  )
+lambda2    <- mxPath(from="F2", to=F2itemNames,values=1, labels = paste("lambda", F2itemNames, sep = "_")  )
+facCor    <- mxPath(from="F1", to= "F2", arrows=2,free=T, values=.3, labels = c("corrs"))
+
+snpMu     <- mxPath(from = "one", to = "snp" , labels = "snpMean")
+snpBeta   <- mxPath(from = "snp", to = latents, labels = paste0("snp", 2, latents), values = 0, free = T)
+snpres    <- mxPath(from = "snp", arrows=2, values=1, free = T, labels = paste("snp", "res", sep = "_"))
+
+resid     <- mxPath(from = c(F1itemNames,F2itemNames), arrows=2, values=1, free = c(fac==0), labels = paste(c(itemNames), "res", sep = "_"))
+facRes    <- mxPath(from=latents, arrows=2,free=F, values=1.0, labels = "facRes")
+covMean   <- mxPath(from = "one", to = covariates, free=FALSE, labels = paste0('data.',covariates)) 
+cov2item  <- mxPath(from = covariates, to = c(c(F1itemNames,F2itemNames)), connect = "all.pairs", labels = paste(rep(covariates, each = length(c(F1itemNames,F2itemNames))), c(F1itemNames,F2itemNames), sep = "_2_"))
+itemMean  <- mxPath(from = 'one', to = c(F1itemNames,F2itemNames), free= c(fac==0), values = 0, labels = paste0(itemNames, "Mean"))
+
+if(maxThr>0) thresh    <- mxThreshold(c(F1itemNames, F2itemNames)[c(fac==1)], nThresh=c(thr[fac==1]), free = T , labels = paste(rep(c(F1itemNames,F2itemNames)[c(fac==1)], each = maxThr), "_Thr_", 1:maxThr, sep = ""), values=mxNormalQuantiles(1))
+
+
+dat       <- mxData(observed=phenoData, type="raw")                                              ## options for min variance/MAF
+
+if(maxThr==0 & fitfun == "WLS") fun      <- mxFitFunctionWLS(allContinuousMethod= "marginals")
+if(maxThr>0 & fitfun == "WLS") fun       <- mxFitFunctionWLS()
+if(fitfun == "FIML") fun                 <- mxFitFunctionML()
+
+
+twoFacPre <- mxModel("TwoFac", type="RAM",
+      manifestVars = c("snp", F1itemNames, F2itemNames),
+      latentVars = c(latents, covariates),
+ 	  lambda1, lambda2, facCor, snpMu, snpBeta, snpres, 
+	  resid, facRes, covMean, cov2item,
+ 	  itemMean, dat, fun  )
+
+
+if(maxThr>0) twoFacPre <- mxModel(twoFacPre, name = "TwoFac", thresh  )
+
+ if(snpFileType == "bgen"){
+ Compute <- mxComputeLoop(list(
+ 	 mxComputeSetOriginalStarts(),
+      mxComputeLoadData("TwoFac", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+ 	 mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+                                ), i=1:nSNP)}
+
+ if(snpFileType == "pgen"){
+ Compute <- mxComputeLoop(list(
+      mxComputeSetOriginalStarts(),
+      mxComputeLoadData("TwoFac", column='snp', path=paste(snpData, snpFileType, sep = "."), method=snpFileType),
+      mxComputeLoadContext(path=paste(snpData, "pvar", sep = "."),column=1:3, sep='\t'),
+      mxComputeTryCatch(mxComputeSequence(list(mxComputeSetOriginalStarts(), mxComputeGradientDescent(),mxComputeStandardError() )  )  ),
+      mxComputeCheckpoint(path=paste(out, "log", sep = "."), standardErrors = TRUE)  #####
+                                ), i=1:nSNP)  }
+
+
+ twoFac <- mxModel(twoFacPre, name = "TwoFac", Compute  )
+
+ twoFacFit <- mxRun(twoFac)
+ summary(twoFacFit)
+
+ }
+
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+##############################################################################################################################################################
+

R/package.R

@@ -0,0 +1,13 @@
+#' Genome-wide Structural Equation Modeling
+#'
+#' @docType package
+#' @name gwsem-package
+#' @aliases gwsem
+#'
+#' @description
+#' Need to write description.
+#'
+#' @useDynLib gwsem, .registration = TRUE
+#' @import methods
+#' 
+NULL

README.md

@@ -0,0 +1,25 @@
+# gwsem
+
+<!-- badges: start -->
+[![Travis build status](https://travis-ci.org/jpritikin/gwsem.svg?branch=master)](https://travis-ci.org/jpritikin/gwsem)
+<!-- badges: end -->
+
+The goal of gwsem is to ...
+
+## Installation
+
+You can install the released version of gwsem from [CRAN](https://CRAN.R-project.org) with:
+
+``` r
+install.packages("gwsem")
+```
+
+## Example
+
+This is a basic example which shows you how to solve a common problem:
+
+``` r
+library(gwsem)
+## basic example code
+```
+

gwsem.Rproj

@@ -18,3 +18,4 @@ StripTrailingWhitespace: Yes
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

src/.gitignore

@@ -0,0 +1,2 @@
+*.o
+*.so