# Multiple-indicator DOC model
# MadhurBain Singh, Nathan Gillespie, Hermine Maes
# March 8, 2024

# TUTORIAL SOLUTIONS  ##########################################################

## ACE Model with the Two Latent True Scores ======================

# ACE - correated A, C and E risks on commmon pathways, no causation 
ACE_fit 	<- mxRun( ACE_model )
summary( ACE_fit )

## Look at the standardized ACE estimates of all six observed variables
var_comp <- cbind(
  round(diag2vec(ACE_fit$algebras$VC$result[1:6,19:24]),2),
  round(diag2vec(ACE_fit$algebras$VC$result[1:6,25:30]),2),
  round(diag2vec(ACE_fit$algebras$VC$result[1:6,31:36]),2))
rownames(var_comp) <- vars 
colnames(var_comp) <- c("A","C","E") 	
var_comp


# 1. Find the best-fitting ACE model (no causation) ============================

## "Global" tests of A, C, and E
## For both latent true scores, as well as all observed indicator variables

### Global AE - correlated background A & E only, no causation =================

## Drop the shared environmental (C) variance components of the latent true scores
## Drop the shared environmental (C) variance components of the observed indicators

AE 				<- mxModel( ACE_fit,name="AE")
AE 				<- omxSetParameters( AE, label=c("Cf_x","covCf","Cf_y"),free=F,values=0)   
AE 				<- omxSetParameters( AE, label=Cs_labs,free=F,values=0)  		
AE_fit  	<- mxRun( AE )
summary( AE_fit )


### Global CE - correlated background C & E only, no causation =================

## Drop the additive genetic (A) variance components of the latent true scores
## Drop the additive genetic (A) variance components of the observed indicators

CE 				<- mxModel( ACE_fit,name="CE") 
CE 				<- omxSetParameters( CE, label=c("Af_x","covAf","Af_y"),free=F,values=0)   
CE 				<- omxSetParameters( CE, label=As_labs,free=F,values=0)  		
CE_fit	  <- mxRun( CE )
summary( CE_fit )

### Global E - correlated background E only, no causation ======================

## Drop the additive genetic (A) variance components of the latent true scores
## Drop the additive genetic (A) variance components of the observed indicators
## Drop the shared environmental (C) variance components of the latent true scores
## Drop the shared environmental (C) variance components of the observed indicators

E 					<- mxModel( AE_fit,name="E") 
E 	  			<- omxSetParameters( E, label=c("Af_x","covAf","Af_y"),free=F,values=0)   
E 				  <- omxSetParameters( E, label=As_labs,free=F,values=0)  		
E_fit     	<- mxRun( E )

summary( E_fit )
## Note the missing Std Errors of EF_x and Ef_y.
## In this model, Ef_x = total variance of Fx
## And we have fixed the total variance of Fx at 1.
## So, Ef_x also, indirectly, becomes fixed at 1; hence, no S.E.

##### Compare ACE with global AE, CE and E covariance models (no causation) ==========

subs 				<- c( AE_fit, CE_fit, E_fit )
comps				<- mxCompare( ACE_fit,subs )
comps
#   base comparison ep minus2LL    df      AIC     diffLL diffdf            p
# 1  ACE       <NA> 39 65628.93 47967 65706.93         NA     NA           NA
# 2  ACE         AE 30 65708.05 47976 65768.05   79.11333      9 2.423902e-13
# 3  ACE         CE 30 65947.37 47976 66007.37  318.43878      9 3.182960e-63
# 4  ACE          E 21 68932.43 47985 68974.43 3303.49753     18 0.000000e+00



## ACE is the best-fitting model
## This models assumes no causation, and so is our "null" model for further causal tests.

nullACE_fit <- ACE_fit

## Now, examine the ACE estimates of X and Y.
## Refit the model with intervals=TRUE
nullACE_fit <- mxRun(nullACE_fit, intervals = T)
round(nullACE_fit$output$confidenceIntervals,3)

#       lbound estimate ubound
# Af_x   0.355    0.459  0.563
# Cf_x  -0.064    0.026  0.113
# Ef_x   0.483    0.515  0.549

# Af_y   0.167    0.255  0.345
# Cf_y   0.219    0.296  0.370
# Ef_y   0.420    0.449  0.480

# covAf  0.093    0.165  0.238
# covCf -0.082   -0.020  0.041
# covEf  0.151    0.173  0.196

## X and Y have very different ACE estimates.
## Afx = 0.46, Cfx = 0.02, Efx = 0.52
## Afy = 0.25, Cfy = 0.30, Efy = 0.45

## Thus, we may be able to fit the DoC model to this pair of phenotypes.

## NB: If the two phenotypes have similar ACE estimates, 
## this pair of phenotypes may not be suitable for DoC modeling!

## Causal inference in the DoC model depends on the different model-expected 
## cross-twin cross-trait covariances, given different directions of causation.

## If the two phenotypes have identical ACE estimates, the model-expected 
## cross-twin cross-trait covariances will end up being very similar, 
## even under different directions of causation.



# 2. Test the DoC model (no background covariance/confounding) ================

## We'll now test the DoC model, i.e., the causal paths between Fx and Fy


### Causation from Fx to Fy. No background ACE covariance ===============

Fx_to_Fy 				<- mxModel( nullACE_fit, name="ACE_byx")
Fx_to_Fy 				<- omxSetParameters( Fx_to_Fy, label=c("covAf","covCf","covEf"),free=FALSE,values=0)  
Fx_to_Fy 				<- omxSetParameters( Fx_to_Fy, label="byx",free=TRUE,values=0.2)  
Fx_to_Fy_fit  	<- mxRun( Fx_to_Fy )
summary( Fx_to_Fy_fit )

### Causation from Fy to Fx. No background ACE covariance ===============

Fy_to_Fx 				<- mxModel( nullACE_fit, name="ACE_bxy")
Fy_to_Fx 				<- omxSetParameters( Fy_to_Fx, label=c("covAf","covCf","covEf"),free=FALSE,values=0)  
Fy_to_Fx 				<- omxSetParameters( Fy_to_Fx, label="bxy",free=TRUE,values=0.2)  
Fy_to_Fx_fit  	<- mxTryHard( Fy_to_Fx )
summary( Fy_to_Fx_fit )


### Bidirectional Causation between Fy to Fx. No background ACE covariance ===============

recip_doc 					<- mxModel( nullACE_fit,name="ACE_recip")
recip_doc 					<- omxSetParameters( recip_doc, label=c("covAf","covCf","covEf"),free=FALSE,values=0)  
recip_doc 					<- omxSetParameters( recip_doc, label=c("byx"),free=TRUE,values=0.2)
recip_doc 					<- omxSetParameters( recip_doc, label=c("bxy"),free=TRUE,values=0.1) 
recip_doc_fit	    	<- mxRun( recip_doc )
summary( recip_doc_fit	)

# Compare the fit of the causal models against the null
subs 						<- c( Fx_to_Fy_fit, Fy_to_Fx_fit, recip_doc_fit )
comps						<- mxCompare( nullACE_fit, subs )
comps
#   base comparison ep minus2LL    df      AIC      diffLL diffdf            p
# 1  ACE       <NA> 39 65628.93 47967 65706.93          NA     NA           NA
# 2  ACE    ACE_byx 37 65631.23 47969 65705.23  2.29420360      2 3.175558e-01
# 3  ACE    ACE_bxy 37 65665.99 47969 65739.99 37.05590365      2 8.982821e-09
# 4  ACE  ACE_recip 38 65628.95 47968 65704.95  0.01610365      1 8.990194e-01

## Reciprocal DOC model model has the best model fit (lowest AIC)*
## *Note that its AIC is only marginally lower than that of Fx_to_Fy DoC model!
## the Fx_to_Fy DOC model fit almost equally well, with one fewer parameter.

## You may select either model to report the causal estimates.
## Obtain the CIs for the freely estimated sources of covariance 
## (causation or background confounding) in the best-fitting model.

best_fit <- recip_doc_fit
best_fit <- mxRun( best_fit, intervals = TRUE)

round(best_fit$output$confidenceIntervals, 3)
##      lbound estimate ubound
# bxy   -0.260   -0.108  0.030
# byx    0.299    0.412  0.530 

## Note that bxy (Y --> X) is not significant.
## So, the causation is likely unidirectional X --> Y

## Can we drop bxy?
## We may directly compare the Fx_to_Fy DoC model with the reciprocal DoC 
## (given non-significant estimates of bxy in the reciprocal model)
mxCompare(recip_doc_fit, Fx_to_Fy_fit)
#        base comparison ep minus2LL    df      AIC diffLL diffdf         p
# 1 ACE_recip       <NA> 38 65628.95 47968 65704.95     NA     NA        NA
# 2 ACE_recip    ACE_byx 37 65631.23 47969 65705.23 2.2781      1 0.1312125

## Thus, the Fx_to_Fy DoC model is the most parsimonious model.
## Even though it has marginally higher AIC than the reciprocal model, 
## the LRT is not significant.

## We should also report the causal estimate from this unidirectional model
## You might have selected the Fx_to_Fy DOC model as the best-ftting one to begin with.
best_fit <- Fx_to_Fy_fit
best_fit <- mxRun( best_fit, intervals = TRUE)
summary( best_fit	)
#               lbound   estimate    ubound 
# byx       0.30533385 0.32526903 0.3449828

## Regardless of which model you keep (unidirectional or reciprocal DoC),
## the DoC model suggests unidirectional causation X --> Y



# BONUS: Test a hybrid model (causation + confounding) =============================

## Notably, DoC assumes no background confounding
## So, in all of the models so far,
## the observed covariance between X and Y (or their indicators) is either 
## due to background confounding
## or, causation
## Some phenotype pairs may have a mix of both!

## In "bivariate" twin models, we can estimate up to **three** sources of covariance
## between the two phenotypes (here, the two latent true scores)

## For example, freely estimate background A covariance, plus reciprocal causation

## Try hybrid model covA, byx, bxy
ACE_hybrid 					<- mxModel( recip_doc_fit,name="Hybrid_covA_byx_bxy")
ACE_hybrid 					<- omxSetParameters( ACE_hybrid, label=c("covAf"),free=TRUE,values=0.2)  
ACE_hybrid_fit	    <- mxRun( ACE_hybrid )
summary( ACE_hybrid_fit	)
# mxCheckIdentification(ACE_hybrid_fit)

# Compare the model fit
# NB: the causation-only model is a restricted, nested model of the hybrid model
mxCompare( ACE_hybrid_fit, recip_doc_fit ) 
# #                base comparison ep minus2LL    df      AIC     diffLL diffdf         p
# 1 Hybrid_covA_byx_bxy       <NA> 39 65628.93 47967 65706.93         NA     NA        NA
# 2 Hybrid_covA_byx_bxy  ACE_recip 38 65628.95 47968 65704.95 0.01610365      1 0.8990194


## Other hybrid models
## covA, covC, byx 
ACE_hybrid2 					<- mxModel( Fx_to_Fy_fit,name="Hybrid_covA_covE_byx")
ACE_hybrid2 					<- omxSetParameters( ACE_hybrid2, label=c("covAf","covEf"),free=TRUE,values=0.2)  
ACE_hybrid2_fit	      <- mxTryHard( ACE_hybrid2 )
summary( ACE_hybrid2_fit	)
# mxCheckIdentification(ACE_hybrid2_fit)

mxCompare( ACE_hybrid2_fit, Fx_to_Fy_fit ) 

# Not differentiable from the other hybrid model above based on model fit statistics
# They both have the same number of parameters
mxCompare( ACE_hybrid2_fit, ACE_hybrid_fit ) 
  
# END. ===========