#
#  OpenMx for ordinal factor analysis of drug data
#

# Load packages

require(psych)
require(OpenMx)
require(polycor)

# Read in original data

drd2data<-read.table('  ',header=T,na.strings='.')
describe(drd2data)



# Ignore correlated status of sibs (naughty naughty) to build big sample size
sib1<-drd2data[,c(8,9,16:18)]
sib2<-drd2data[,c(10,11,19:21)]
sib3<-drd2data[,c(12,13,22:24)]
sib4<-drd2data[,c(14,15,25:27)]
names(sib1)<-c("SNP1","SNP2","Stimulants","Tranquilizers","Marijuana")
names(sib2)<-names(sib1)
names(sib3)<-names(sib1)
names(sib4)<-names(sib1)

sibsasInd<-rbind(sib1,sib2,sib3,sib4)

# Retain only those cases who are not missing on all variables
sibsWithData<-subset(sibsasInd,!(is.na(Stimulants) & is.na(Tranquilizers) & is.na(Marijuana)))
sibsWithSNPs<-subset(sibsWithData,!is.na(SNP1) & !is.na(SNP2))
# Take a look at the correlations in two ways
cor(sibsWithData,use="pairwise")
tetrachoric(sibsWithData)

# Now try to build and fit factor model
# Step 5: make the ordinal variables into R factors
ordinalData <- cbind(mxFactor(sibsWithSNPs[,3:5],levels=c(0:1)),sibsWithSNPs[,1:2])
nVariables <- 3
nFactors <- 1
nThresholds <- 1
obsVarNames <- names(ordinalData)[1:3]

# Fit Factor Model with Raw Ordinal Data and Matrices Input
# -----------------------------------------------------------------------
oneFactorThresholdModel <- mxModel('oneFactorThresholdModel',

    mxMatrix( type="Full", nrow=nVariables, ncol=nFactors, free=TRUE, values=0.2, lbound=-.99, ubound=.99, name="facLoadings" ),
    mxMatrix( type="Unit", nrow=nVariables, ncol=1, name="vectorofOnes" ),
    mxMatrix( type="Zero", nrow=1, ncol=nVariables, name="Zeroes" ),
    mxMatrix( type="Full", nrow=nThresholds, ncol=nVariables, free=TRUE, values=.2, dimnames=list(c(), obsVarNames),
              name="thresholdDeviations" ),
    mxMatrix( type="Lower", nrow=nThresholds, ncol=nThresholds, free=FALSE, values=1, name="unitLower" ),

    mxAlgebra( expression=vectorofOnes - (diag2vec(facLoadings %*% t(facLoadings))) , name="resVariances" ),
    mxAlgebra( expression=facLoadings %*% t(facLoadings) + vec2diag(resVariances), name="expCovariances" ),
    mxAlgebra( expression=unitLower %*% thresholdDeviations, name="expThresholds" ),
    mxAlgebra( expression=Zeroes, name="expMeans" ),

    mxData( observed=ordinalData, type='raw' ),
    mxFIMLObjective( covariance="expCovariances", means="expMeans", dimnames=obsVarNames, thresholds="expThresholds" )
)

oneFactorThresholdFit <- mxRun(oneFactorThresholdModel, suppressWarnings=TRUE)
summary(oneFactorThresholdFit)


# -----------------------------------------------------------------------
# Fit Factor Model SNPs influencing latent trait
# -----------------------------------------------------------------------

# SNPs influence latent trait
# Delete the old expMeans
oneFactorThresholdModel$expMeans <- NULL

oneFactorThresholdModelSNP <- mxModel(oneFactorThresholdModel,
    mxMatrix( type="Full", nrow=2, ncol=1, free=F, values=0, labels=c("data.SNP1","data.SNP2"), name="SNPs" ),
    mxMatrix( type="Full", nrow=1, ncol=1, free=T, name="Beta"),
    mxAlgebra( ((sum(SNPs) * Beta) %x% t(facLoadings)), name="expMeans")
)

oneFactorThresholdSNPFit <- mxRun(oneFactorThresholdModelSNP, suppressWarnings=TRUE)
summary(oneFactorThresholdSNPFit)

# SNPs influence each substance
# Delete the old expMeans
oneFactorThresholdModel$expMeans <- NULL

oneFactorThresholdModelSNPres <- mxModel(oneFactorThresholdModelSNP,
    mxMatrix( type="Full", nrow=1, ncol=nVariables, free=T, name="Beta"),
    mxAlgebra( sum(SNPs) %x% Beta, name="expMeans")
)

oneFactorThresholdSNPresFit <- mxRun(oneFactorThresholdModelSNPres, suppressWarnings=TRUE)
summary(oneFactorThresholdSNPresFit)


# -----------------------------------------------------------------------

# Fit Latent Class Model with Raw Ordinal Data and Matrices Input
# -----------------------------------------------------------------------

nVar<-3
nClass<-2

# -----------------------------------------------------------------------
# Create an MxModel object
# -----------------------------------------------------------------------
class1 <- mxModel("Class1", 
    mxMatrix( type="Iden", nrow=nVar, ncol=nVar, name="S" ),
    mxMatrix( type="Full", nrow=1, ncol=nVar, values=c(0.1,0.6,0.9), free=T, name="Thresholds" ),
    mxMatrix( type="Zero", nrow=1, ncol=nVar, name="Zeroes" ),
    mxData( observed=ordinalData, type='raw' ),
	mxFIMLObjective(covariance="S", means="Zeroes", dimnames=names(ordinalData[,1:3]), thresholds="Thresholds", vector=T)
)

# Make a class 2 model that looks pretty much the same as class 1 but has a different name
class2 <- mxModel(class1, name="Class2" ) 

# Nudge class2 endorsement probs away from those of class 1
class2@matrices$Thresholds@values[]<-0

# make a matrix of class probabilities
classP <- mxMatrix("Full", nClass, 1, free=c(TRUE, FALSE), values=2, lbound=0.001, labels = c("p1", "p2"), name="Props")

# standardize them
classS <- mxAlgebra(Props%x%(1/sum(Props)), name="classProbs")

# build the algebra to compute -2lnL where L is sum of (wi * L(items | class=i) ) with wi being class i's memb prob 
algObj <- mxAlgebra(-2*sum( log(classProbs[1,1]%x%Class1.objective + classProbs[2,1]%x%Class2.objective)), name="mixtureObj")

# make a mxAlgebraObjective object
obj <- mxAlgebraObjective("mixtureObj")

# bundle all the bits together into one mxModel       
lcaModel <- mxModel("Latent Class Model", 
                    mxData( observed=ordinalData[,1:3], type='raw' ),class1, class2, classP, classS, algObj, obj )

lcaModelFit <- mxRun(lcaModel, suppressWarnings=TRUE)
summary(lcaModelFit)
# Take a look at some of the output
lcaModelFit$classProbs




# Modify model to allow SNPs to affect classProbs
# Trick is to put Definition variables into a matrix to evaluate classProps accordingly

SNPsums <- as.matrix(ordinalData[,4]) + as.matrix(ordinalData[,5])
U <- mxMatrix( "Unit", nrow=dim(ordinalData)[1], ncol=1, free=F, name="U")
U22 <- mxMatrix( "Unit", nrow=2, ncol=2, name="U22")
SNP <- mxMatrix( type="Full", nrow=dim(ordinalData)[1], ncol=1, free=F, values=SNPsums, name="SNP")
Beta <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), name="Beta")
BetaSNP <- mxAlgebra( SNP %x% Beta, name="BetaSNP")
UProp <- mxAlgebra( U %x% t(Props), name="UProp") 

classS <- mxAlgebra( (UProp + BetaSNP) / ((UProp + BetaSNP) %*% U22), name="classProbs")
algObj <- mxAlgebra(-2*sum( log(classProbs[,1] * Class1.objective + classProbs[,2] * Class2.objective)), name="mixtureObj")

lcaModelSNP <- mxModel("Latent Class Model", 
                    mxData( observed=ordinalData[1:3,], type='raw' ),class1, class2, classP, classS, SNP, Beta, BetaSNP, U, U22, UProp, algObj, obj )
lcaModelSNPFit <- mxRun(lcaModelSNP)
summary(lcaModelSNPFit)


### End