# Check identicaion of reciprocal causation model fitted to bivariate data with genetic liabilties # measured directly (i.e. genomically) # Using Open Mx with expected "observed" covariance matrix generated deterministically # ****** Model assumes G1,GC and G2 correlated (a,b,c ne 0), g1,g2,h1,h2 ne 0; bidirectional causation p,q ne 0 # f1 is zero and correlated residuals (r>0) # ***** Part 1 ******** # Generate "observed" standardized covariance matrices (Sigma) under full (e.g. "Mendelian randomization") causal model # Model assumes quantitative genetic liabilities can be measured without error and resolved into outcome specific and # common effects shared by both (easier said than done!) # Supply model parameters (refer to accompanying powerpoint slides) # Initial genetic parameters and residual paths are standardized to unit total phenotypic variance. # After reciproacal causation the total variance is changed (diagonal of S). # Sigma is then rescaled to unit total variance # Parameters: g1 path from genetic effects specific to P1; g2 -do- P1; h1 path from shared genes to P1; # f1 path from genetic effects on P1 to P2 (usually assumed to be zero); h2 -do- P2; # a = correlation of two sets of specific genes (e.g. stratificaion); b -do- G1,Gc; c -do- Gc,G2; # r1, r2 residual paths to P1, P2; r=correlation between residuals; # p and q = reciprocal interactions p=P1->P2; q= -do- P2->P1. b<- 0.2 #r(G1,Gc) c<- 0.2 #r(G2,Gc) a<- 0.1 #r(G1,G2) g1<- 0.4 #p(G1,P1) h1<- 0.3 #p(Gc,P1) f1<- 0.0 #p(G1,P2) g2<- 0.2 #p(G2,P2) h2<- 0.4 #p(Gc,P2) p<- 0.2 #p(P1,P2) q<- 0.3 #p(P2,P1) r<- 0.2 #r(R1,R2) # Residual paths obtained by difference assuming G, R and P are iniially standardised to unit variance (ignoring p and q). r1<- sqrt(1-(g1^2+h1^2+2*b*h1*g1)) r2<- sqrt(1-(f1^2+g2^2+h2^2+2*c*h2*g2+2*b*f1*h2+2*+2*a*f1*g2)) # Order of variables in matrices: G1, Gc, G2, R1, R2, P1, P2 # Generate Initial matrices between measured genetic, residual, and phenotypic variables: G, R and P R<- matrix(nrow=7,ncol=7, data=c(1,b,a,0,0,0,0, b,1,c,0,0,0,0, a,c,1,0,0,0,0, 0,0,0,1,r,0,0, 0,0,0,r,1,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0),byrow=TRUE) # Generate Beta matrix from parameters Beta<- matrix(nrow=7,ncol=7 , data=c(0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, g1,h1,0,r1,0,0,q, f1,h2,g2,0,r2,p,0),byrow=TRUE) I<- matrix(nrow=7,ncol=7, data= c(1,0,0,0,0,0,0, 0,1,0,0,0,0,0, 0,0,1,0,0,0,0, 0,0,0,1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1) ,byrow=TRUE) # Obtain expected "observed" covariance matrix for infinite sample size TT<-solve(I-Beta) # Obtain matrix to remove unwanted rows and columns (4 qnd 5 correspond to latent residuals) T<- matrix(nrow=5,ncol=7,data=c(1,0,0,0,0,0,0, 0,1,0,0,0,0,0, 0,0,1,0,0,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1),byrow=TRUE) # Compute covariance matrix generated by model parameters for infinite sample size S<-matrix(T%*%TT%*%R%*%t(TT)%*%t(T),nrow=5) # Standardize S Delta<-solve(sqrt(diag(diag(S)))) Sigma<-matrix(Delta%*%S%*%Delta,nrow=5, dimnames=list(c("g1","gc","g2","p1","p2"),c("g1","gc","g2","p1","p2"))) # Round Sigma to get rid of asymmetry due to rounding errors Sigma<-signif(Sigma,digits=6) Sigma # ******* Part 2 ********* # Fit model to recover parameter values # Set up Mx model #### IMPORTANT NOTE ABOUT SCALING OF RESULTS FOR COMPARISON WITH SUPPLIED VALUES ##################### # a,b,c and r are correct directly from Mx results. g1,g2,h1,h2,r1 and r2 are rescaled subsequently for comparison # with supplied values. # p and q have to be rescaled to recover values supplied to original Beta matrix ###################################################################################################### require(OpenMx) # Correlation matrix of genetic measures and residuals (phenotypes included with zero variance at this stage) Rmat<-mxMatrix(type="Symm",nrow=7,ncol=7, values=c(1,0,0,0,0,0,0, 1,0,0,0,0,0, 1,0,0,0,0, 1,0.1,0,0, 1,0,0, 0,0, 0), free=c(FALSE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE, FALSE,TRUE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,TRUE,FALSE,FALSE, FALSE,FALSE,FALSE, FALSE,FALSE, FALSE), labels=c(NA,"b","a",NA,NA,NA,NA, NA,"c",NA,NA,NA,NA, NA,NA,NA,NA,NA, NA,"r",NA,NA, NA,NA,NA, NA,NA, NA),name="R") Bmat<-mxMatrix(type="Full",nrow=7,ncol=7, byrow=TRUE, values=c(0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0.,0.2,0,0.93,0,0,0.0, 0.,0.3,0.0,0,0.95,0.0,0), free=c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, TRUE,TRUE, FALSE,TRUE ,FALSE,FALSE,TRUE, FALSE,TRUE, TRUE,FALSE,TRUE,TRUE ,FALSE), labels=c(NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA, "g1","h1",NA,"r1",NA,NA,"q", "f1","h2","g2",NA,"r2","p",NA),name="Beta") Imat<-mxMatrix(type="Iden",nrow=7,ncol=7,name="I") TTmat<-mxAlgebra(expression=solve(I-Beta),name="TT") Tmat<-mxMatrix(type="Full",nrow=5,ncol=7,values=c(1,0,0,0,0,0,0, 0,1,0,0,0,0,0, 0,0,1,0,0,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1),free=FALSE, byrow=TRUE,name="T") Umat<-mxMatrix(type="Unit",nrow=5,ncol=1,name="One") ESigmamat<-mxAlgebra(expression=(T%*%TT)%&%R,name="ESigma") # Extract diagonal of ESigma Dmat<-mxAlgebra(expression=diag2vec(ESigma),name="Delta") # Calculate scale parameters for genetic and residual paths Scale1<-mxAlgebra(expression=sqrt(Beta[6,1]^2+Beta[6,2]^2+2*Beta[6,1]*Beta[6,2]*R[1,2]+Beta[6,4]^2),name="S1") Scale2<-mxAlgebra(expression=sqrt(Beta[7,1]^2+Beta[7,2]^2+Beta[7,3]^2+Beta[7,5]^2+ 2*Beta[7,1]*Beta[7,2]*R[1,2]+2*Beta[7,1]*Beta[7,3]*R[1,3]+2*Beta[7,2]*Beta[7,3]*R[2,3]),name="S2") # Rescale paths to unit total initial variance (to check paths are same as raw paths in simulation) g1<-mxAlgebra(expression=Beta[6,1]/S1,name="G1") h1<-mxAlgebra(expression=Beta[6,2]/S1,name="H1") f1<-mxAlgebra(expression=Beta[7,1]/S2,name="F1") g2<-mxAlgebra(expression=Beta[7,3]/S2,name="G2") h2<-mxAlgebra(expression=Beta[7,2]/S2,name="H2") r1<-mxAlgebra(expression=Beta[6,4]/S1,name="R1") r2<-mxAlgebra(expression=Beta[7,5]/S2,name="R2") # Constraint expected variances to unity Con1<-mxConstraint(expression=Delta==One,name="C1") #Check algebra (not needed when working) #Model<-mxModel(name="MyAlgebra", Rmat,Bmat,Imat,TTmat,Tmat,Dmat,Umat,ESigmamat,Con1) #Eval<-mxRun(Model) #Eval@algebras # Supply data, objective function and fit model data<-mxData(observed=Sigma,type="cov",numObs=1000) objective=mxMLObjective(covariance="ESigma",dimnames=c("g1","gc","g2","p1","p2")) MendelModel<-mxModel(name="MendelianRandomization",Rmat,Bmat,Imat,TTmat,Tmat,ESigmamat,Umat,Dmat, Scale1,Scale2,g1,h1,f1,r1,g2,h2,r2,Con1,data,objective) MendelFit<-mxRun(MendelModel) summary(MendelFit) MendelFit@algebras