However, I got a warning like this: Is this something important? Finally, I am not sure if I have set the setting values correctly. How can I calculate?
Warning messages:
1: In matrix(labels, nrow, ncol, byrow = byrow) :
data length [3] is not a sub-multiple or multiple of the number of columns [2]
2: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'values' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh("th", vars, nth), name = "thinG")
3: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh("th", vars, nth), name = "thinG")
It looks as though you're trying to put something of length 3 into a 1x2 matrix. Try inspecting thingG by printing it in the R console, and see if its 'labels' matrix looks correct. If I understand you correctly, you have two continuous traits and one ordinal trait, so you have a total of two ordinal variables, right?
Yes that's true. Now I checked the label part as you said, it appears as 1x2 matrix. How can I fix this? I have 2 continuous trait (cape and no_2_FA) and one dummy coded variable (cannabis_dummy 0:non-user 1:user).
Thank you problem solved. I tried setting up a trivariate cholesky model using the same variables. I also tried adding tobacco and alcohol to covariate, but I am getting an error like this. I think I'm making a mistake while creating the matrix. How can I fix this. My script is as follows:
library(OpenMx)library(psych);library(polycor)source("miFunctions.R")library(readxl)
TwinData<-read_excel("/Users/handeezgiatmaca/Desktop/cannabis_dataset_wide.xlsx")#read the excel filehead(TwinData)
describe(TwinData)names(TwinData)#Regress out age + sex from continuous data (DWI)#Twin 1residuals(lm(TwinData$no_2_FA1 ~ TwinData$age1+ TwinData$gender_dummy1, na.action="na.exclude"))-> TwinData$r_no_2_FA1
#Twin 2residuals(lm(TwinData$no_2_FA2 ~ TwinData$age2+ TwinData$gender_dummy2, na.action="na.exclude"))-> TwinData$r_no_2_FA2
#Regress out age + sex from continuous data (Cape)#Twin 1residuals(lm(TwinData$cape1 ~ TwinData$age1+ TwinData$gender_dummy1, na.action="na.exclude"))-> TwinData$r_cape1
#Twin 2residuals(lm(TwinData$cape2 ~ TwinData$age2+ TwinData$gender_dummy2, na.action="na.exclude"))-> TwinData$r_cape2
#Regress out age + sex from continuous data (alcohol)#Twin 1residuals(lm(TwinData$alcohol_week1 ~ TwinData$age1+ TwinData$gender_dummy1, na.action="na.exclude"))-> TwinData$r_alcohol_week1
#Twin 2residuals(lm(TwinData$alcohol_week2 ~ TwinData$age2+ TwinData$gender_dummy2, na.action="na.exclude"))-> TwinData$r_alcohol_week2
#Regress out age + sex from continuous data (tobacco)#Twin 1residuals(lm(TwinData$tobacco_daily1~ TwinData$age1+ TwinData$gender_dummy1, na.action="na.exclude"))-> TwinData$r_tobacco_daily1
#Twin 2residuals(lm(TwinData$tobacco_daily1 ~ TwinData$age2+ TwinData$gender_dummy2, na.action="na.exclude"))-> TwinData$r_tobacco_daily2
names(TwinData)#Check SD for Twin 1 & Twin2sd(TwinData$r_no_2_FA1,na.rm=T)sd(TwinData$r_no_2_FA2,na.rm=T)sd(TwinData$r_cape1,na.rm=T)sd(TwinData$r_cape2,na.rm=T)sd(TwinData$r_alcohol_week1,na.rm=T)sd(TwinData$r_alcohol_week2,na.rm=T)sd(TwinData$r_tobacco_daily1,na.rm=T)sd(TwinData$r_tobacco_daily2,na.rm=T)# Multiply to increase SD
TwinData$r_no_2_FA1<-(TwinData$r_no_2_FA1 *24)sd(TwinData$r_no_2_FA1,na.rm=T)
TwinData$r_no_2_FA2 <-(TwinData$r_no_2_FA2*16.7)sd(TwinData$r_no_2_FA2,na.rm=T)
TwinData$r_cape1 <-(TwinData$r_cape1/14.7)sd(TwinData$r_cape1,na.rm=T)
TwinData$r_cape2 <-(TwinData$r_cape2/16.8)sd(TwinData$r_cape2,na.rm=T)
TwinData$r_alcohol_week1 <-(TwinData$r_alcohol_week1/2.3)sd(TwinData$r_alcohol_week1,na.rm=T)
TwinData$r_alcohol_week2 <-(TwinData$r_alcohol_week2/2.6)sd(TwinData$r_alcohol_week2,na.rm=T)
TwinData$r_tobacco_daily1 <-(TwinData$r_tobacco_daily1/6.5)sd(TwinData$r_tobacco_daily1 ,na.rm=T)
TwinData$r_tobacco_daily2 <-(TwinData$r_tobacco_daily2/6.9)sd(TwinData$r_tobacco_daily2,na.rm=T)#Check Normality hist(TwinData$r_no_2_FA1)hist(TwinData$r_no_2_FA2)hist(TwinData$r_cape1)hist(TwinData$r_cape2)# Select Continuous Variables
varsc <-c('r_no_2_FA','r_Cape')# list of continuous variables names
nvc <-2# number of continuous variables
ntvc <- nvc*2# number of total continuous variables
conVars <-paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="")# Select Ordinal Variables
nth <-1# number of thresholds
varso <-c('cannabis_dummy')# list of ordinal variables names
nvo <-1# number of ordinal variables
ntvo <- nvo*2# number of total ordinal variables
ordVars <-paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData <- TwinData
#for (i in ordVars) { ordData[,i] <- cut(twinData[,i], #breaks=quantile(ordData[,i],(0:(nth+1))/(nth+1),na.rm=TRUE), labels=c(0:nth)) }# Select Variables for Analysis
vars <-c('r_no_2_FA','r_cape','cannabis_dummy')# list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2# number of total variables
oc <-c(0,0,1)# 0 for continuous, 1 for ordinal variables
selVars <-paste(vars,c(rep(1,nv),rep(2,nv)),sep="")# Select Covariates for Analysis
covVars <-c('r_tobacco_daily1', "r_tobacco_daily2" , "r_alcohol_week1","r_alcohol_week2")
nc <-4# number of covariates # number of covariates# Select Data for Analysis
mzData <-subset(TwinData, zygosity=="MZ", c(selVars, covVars))
dzData <-subset(TwinData, zygosity=="DZ", c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$cannabis_dummy1 <- mxFactor(mzDataF$cannabis_dummy1, levels=c(0,1))
mzDataF$cannabis_dummy2 <- mxFactor(mzDataF$cannabis_dummy2, levels=c(0,1))
dzDataF$cannabis_dummy1 <- mxFactor(dzDataF$cannabis_dummy2, levels=c(0,1))
dzDataF$cannabis_dummy2 <- mxFactor(dzDataF$cannabis_dummy2, levels=c(0,1))# Generate Descriptive Statisticsapply(mzData,2,myMean)apply(dzData,2,myMean)# Set Starting Values
frMV <-c(TRUE,FALSE)# free status for variables
frCvD <-diag(frMV,ntv,ntv)# lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)]<- TRUE # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)]<- TRUE # lower bounds for above diagonal elements
frCv <-matrix(as.logical(frCvD),4)
svMe <-c(0,0,0)# start value for means
svPa <- .1 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv))# start values for diagonal of covariance matrix
svPe <- .1 # start value for path coefficient for e
svPeD <- vech(diag(svPe,nv,nv))# start values for diagonal of covariance matrix
lbPa <- .0001 # start value for lower bounds
lbPaD <-diag(lbPa,nv,nv)# lower bounds for diagonal of covariance matrix
lbPaD[lower.tri(lbPaD)]<--10# lower bounds for below diagonal elements
lbPaD[upper.tri(lbPaD)]<- NA # lower bounds for above diagonal elements
svLTh <--1# start value for first threshold
svITh <-1# start value for increments
svTh <-matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv)# start value for thresholds
lbTh <-matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv)# lower bounds for thresholds# Create Labels
labMe <-paste("mean",vars,sep="_")
labTh <-paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_")
labBe <- labFull("beta",nc,1)# ------------------------------------------------------------------------------# PREPARE MODEL# ACE Model# Create Matrices for Covariates and linear Regression Coefficients
deftobacco <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.r_tobacco_daily1", "data.r_tobacco_daily1","data.r_tobacco_daily1","data.r_tobacco_daily2","data.r_tobacco_daily2", "data.r_tobacco_daily2"), name="tobacco")
pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=6, free=TRUE, values=.01, label=c("b11","b12","b13","b11","b12","b13"), name="b1")
defalcohol <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.r_alcohol_week1", "data.r_alcohol_week1","data.r_alcohol_week1","data.r_alcohol_week2","data.r_alcohol_week2", "data.r_alcohol_week2"), name="alcohol")
pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23","b21", "b22", "b23"), name="b2")# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMe, name="meanG")
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG")
expMean <- mxAlgebra(expression= meanG +cbind(t(b1%*%tobacco),t(b1%*%tobacco))+cbind(t(b2%*%alcohol),t(b2%*%alcohol)) , name="expMean")
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinG")
inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc")
threG <- mxAlgebra(expression= inc %*% thinG, name="threG")# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a")
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c")
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e")# Create Algebra for Variance Comptwonts
covA <- mxAlgebra(expression=a %*%t(a), name="A")
covC <- mxAlgebra(expression=c%*%t(c), name="C")
covE <- mxAlgebra(expression=e %*%t(e), name="E")# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra(expression= A+C+E, name="V")
covMZ <- mxAlgebra(expression= A+C, name="cMZ")
covDZ <- mxAlgebra(expression=0.5%x%A+C, name="cDZ")
expCovMZ <- mxAlgebra(expression=rbind(cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ")
expCovDZ <- mxAlgebra(expression=rbind(cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ")# Create Algebra for Standardization
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra(expression=solve(sqrt(I*V)), name="iSD")# Calculate genetic and environmental correlations
corA <- mxAlgebra(expression=solve(sqrt(I*A))%&%A, name ="rA")#cov2cor()
corC <- mxAlgebra(expression=solve(sqrt(I*C))%&%C, name ="rC")
corE <- mxAlgebra(expression=solve(sqrt(I*E))%&%E, name ="rE")# Constrain Variance of Binary Variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1")
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc")
var1 <- mxConstraint(expression=diag2vec(Oc%&%V)==Unv1, name="Var1")# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzDataF, type="raw")
dataDZ <- mxData( observed=dzDataF, type="raw")# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
funML <- mxFitFunctionML()# Create Model Objects for Multiple Groups
pars <-list(pathB1,pathB2, meanG, thinG, matI, invSD, matUnv, matOc,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
defs <-list(deftobacco,defalcohol)
modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup(c("MZ","DZ"))# Create Algebra for Variance Components
rowVC <-rep('VC',nv)
colVC <-rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC <- mxAlgebra(expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))# Create Confidence Interval Objects
ciACE <- mxCI("VC[1,1]")# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" )# Build Model with Confidence Intervals
modelACE <- mxModel("twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE )# ------------------------------------------------------------------------------# RUN MODEL# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T)
sumACE <-summary( fitACE )# Compare with Saturated Model
mxCompare( fit, fitACE )# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEsts(fitACE)# ------------------------------------------------------------------------------# RUN SUBMODELS# Run AE model
modelAE <- mxModel( fitACE, name="twoAEja")
modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0)
fitAE <- mxRun( modelAE, intervals=T)
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)# Run CE model
modelCE <- mxModel( fitACE, name="twoCEja")
modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0)
fitCE <- mxRun( modelCE, intervals=T)
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
fitEsts(fitCE)# Run E model
modelE <- mxModel( fitAE, name="twoEja")
modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0)
fitE <- mxRun( modelE, intervals=T)
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)# Print Comparative Fit Statistics
mxCompare( fitACE, nested <-list(fitAE, fitCE, fitE))
The error is like this:
># Run ACE Model> fitACE <- mxRun( modelACE, intervals=T)
Error in value[rows[[i]], cols[[i]]]<- startValue :
incorrect number of subscripts on matrix
Try running traceback() after you get the error, and report what traceback() says.
Also, why are you adjusting for age and sex outside your MxModel? You can use age and sex as definition variables, just like tobacco and alcohol usage.
Printer-friendly version
It looks as though you're trying to put something of length 3 into a 1x2 matrix. Try inspecting
thingGby printing it in the R console, and see if its 'labels' matrix looks correct. If I understand you correctly, you have two continuous traits and one ordinal trait, so you have a total of two ordinal variables, right?Yes that's true. Now I checked the label part as you said, it appears as 1x2 matrix. How can I fix this? I have 2 continuous trait (cape and no_2_FA) and one dummy coded variable (cannabis_dummy 0:non-user 1:user).
Thank you for your help!
Try defining
svTh,lbTh, andthinGdifferently, as follows:Thank you problem solved. I tried setting up a trivariate cholesky model using the same variables. I also tried adding tobacco and alcohol to covariate, but I am getting an error like this. I think I'm making a mistake while creating the matrix. How can I fix this. My script is as follows:
The error is like this:
Thank you!
Try running
traceback()after you get the error, and report whattraceback()says.Also, why are you adjusting for age and sex outside your MxModel? You can use age and sex as definition variables, just like tobacco and alcohol usage.