来源网站：http://dbtemp.blogspot.com/2011/08/sem-approach-to-twin-data.html
对于行为遗传学分析来说，我们希望能够衡量一个比例，即人与人之间性状的相似性中有多少可以用基因的相似性来解释。在最简单的双生子分析中，我们的观测变量包括两个协方差矩阵：所有MZ同卵双胞胎一个、所有DZ异卵双胞胎一个。
针对某一性状（变量）而言，协方差矩阵可以表示双胞胎对之间该性状的协方差（可以理解为非标准化的相关系数，反应有单位的共变程度）。
这个简单模型中最重要的假设是：双胞胎之间的形状的相似性（差异）可以用遗传、共有环境（包括出生前胎内环境、家庭环境、社会经济地位、住宅区；使得双胞胎彼此更像）和他们特有的环境（这些环境因素不共有，包括一些特殊的经历[疾病、受伤等]、随机的生物效应、自身对感知的差别、以及测量的误差；造成双胞胎之间的差异）来解释。
一个容易混淆的点是，当我们谈论个体之间遗传相似性的时候，我们认为：同卵双胞胎MZ具有完全相同遗传物质，而异卵双胞胎DZ之间有50%的相同遗传物质。这其实很难理解的（确实），因为我们也常常听到有人说我们人类和大猩猩基因相似性高达98%。这里不如这么理解：在行为遗传学分析中，我们只关注那些在不同人上有不同等位基因形式的基因，这些基因可能可以解释我看观察得到的个体差异。这些具有多态性的基因仅仅占人类整个基因组中很小的一部分。
如果我们暂时只关心加性基因的效应（additive genetic influence），则对于异卵双胞胎可以表示如下图所示的标准ACE模型。
如图所示，每个潜变量都有相对应的路径，我们的目标便是估计出这些路径。这里需要注意的是，twin1和twin2，从潜变量发出的路径我们用了相同的字母，这说明在这个模型中，我们假设对于两个双胞胎来说，基因/环境的影响作用的大小是没有差异的（？译者理解：潜变量本身是有差别的，但是潜变量对与性状的作用是相同的）。这样，也就使得twin1和twin2之间是没有顺序的（是随机的），this is a valid assumption。
上图中，两个潜变量C之间的路径设置为1，即，twin1和twin2从定义上说环境因素是一样的；而两个潜变量A之间的路径被设置为了0.5，也就是，平均来说，DZ双胞胎之间50%的等位基因是相同。而潜变量E之间没有任何设定关系。
- 使用之前的路径追踪准则对模型进行分析，可以估计出：
DZ模型中twin1和twin2的协方差：
CovDZ = 0.5*a2 + 1*c2
总变异方差：
VarDZ = a2 + c2+ e2
- MZ模型与DZ模型很类似，仅仅区别在潜变量A之间的路径系数为1：
模型中twin1和twin2的协方差：
CovMZ = 1*a2 + 1*c2
总变异方差与DZ模型的相同
如果将变量进行标准化，则对于每个双胞胎来说，总的变异为1。则标准化以后估计得到的a2就可以作为加性基因能解释的变异所占的百分比，也就是遗传度。（敲黑板！划重点啊）
上述公式告诉我们，我们可以通过DZ和MZ双胞胎对之间twin1和twin2的相关性对a2 、 c2和e2的大小进行估计。所以我们知道twin1和twin2之间的相关，则可以直接得到遗传度：
因此，最简单的算术方法即为：
(CorMZ-CorDZ) = 0.5 * a2
从而推出：
a2 = 2 *(CorMZ-CorDZ)
c2 = CorMZ - a2
e2 = 1 - a2 - c2
既然这么简单就能计算遗传度，那为什么我们整了这么复杂的模型拟合呢？这里有三个原因。
1.模型拟合是更严格的检验
模型拟合对于违反模型假设的情况是十分敏感的。通常，模型拟合估计的是协方差而不是相关系数，因此我们可以确定，MZ和DZ双胞胎对之间的变异是可比较的（译者：协方差是带单位的，相关不带）。模型假设变异是相同的（VarDZ=VarMZ），所以当不满足假设时，模型便会拟合的不好。
2.模型拟合不仅可以估计路径
还可以同时估计出参数的置信区间，以及可以比较嵌套模型之间的拟合程度以获得额外的信息。
3.可以更容易扩展模型使其适用于多变量分析、或者是更加复杂的模型
在OpenMx中对双生子数据进行建模
在OpenMx中有两种建模的方法：
1.矩阵代数表示
2.RAM路经方法
两种方法都是可以的，但是omxGraphviz的指令只能针对RAM路经法进行绘图。
对双生子建模的老方法一般用MZ和DZ的协方差矩阵作为模型输入，但是现在一般也直接利用被试的原始数据（生数据raw data）当作输入，这样一方面可以更方便使用者，另外一方面我们就可以把确实值包含在内以扩大样本量。事实上，OpenMx的FIML优化方法中，我们除了估计协方差矩阵还会对均值进行估计。
数据生成
下面是一个例子🌰，用了一组模拟的1000对双生子数据，预计a2= 0.5；c2=0.3; e2=0.2
#----------------------------------------------------------------------------------------------------------------------
# Twin simulate
# DVM Bishop, 11th March 2010, Based on script in OpenMXUserGuide, p 15
#----------------------------------------------------------------------------------------------------------------------
require(OpenMx) # 这个是必须的，加载OpenMx包
#------------------下面是模拟数据需要的
require(MASS) # needed for multivariate random number generation
set.seed(200) # specified seed ensures same random number set generated on each run
#-----------设定a、c、e的值
mya2<-0.5 #Additive genetic variance component (a squared)
myc2<-0.3 #Common environment variance component (c squared)
mye2<-1-mya2-myc2 #Specific environment variance component (e squared)
#-----------给定要生成数据的相关系数
my_rMZ <-mya2+myc2 # correlation between MZ twin1 and twin2
my_rDZ <- .5*mya2+myc2 # correlation between DZ twin 1 and twin 2
#-----------利用mvrnorm函数生成数据，MZ、DZ共1000对数据，均值均为0，
#-----------被试之间的协方差矩阵为CovMZ=[1 0.8; 0.8 1];CovDZ=[1 0.55; 0.55 1]
myDataMZ <- mvrnorm (1000, c(0,0), matrix(c(1,my_rMZ,my_rMZ,1),2,2))
myDataDZ <- mvrnorm (1000, c(0,0), matrix(c(1,my_rDZ,my_rDZ,1),2,2))
#-----------给刚刚生成的双生子数据加上标签，并且给出数据的基本统计信息
colnames(myDataMZ) <- c('twin1', 'twin2') # assign column names
colnames(myDataDZ) <- c('twin1', 'twin2')
summary(myDataMZ)
summary(myDataDZ)
colMeans(myDataMZ,na.rm=TRUE)
#----------na.rm 表示忽略 NA值 (非数值)
colMeans(myDataDZ,na.rm=TRUE)
cov(myDataMZ,use="complete")
#----------"complete" 表示使用所有数据，没有确实值
cov(myDataDZ,use="complete")
#---------- 对 MZ 和 DZ数据画散点图
split.screen(c(1,2)) # split display into two screens side by side
# (use c(2,1) for screens one above the other)
screen(1)
plot(myDataMZ,main='MZ') # main specifies overall plot title
screen(2)
plot(myDataDZ, main='DZ')
#use drag and drop to resize the plot window if necessary
#----------把所有数据串起来作为最后生成的数据
alltwin=cbind(myDataMZ,myDataDZ)
colnames(alltwin)=c("MZ_twin1","MZ_twin2","DZ_twin1","DZ_twin2")
write.table(alltwin,"mytwinfile")
# 保存到R当前的路经下 "mytwinfile"
#--------------------------------------------------------------------------------------------
在有缺失值的情况下，上述代码可以做以下改变：
nudat=edit(myDataDZ) #拷贝变量myDataDZ并指定那些值是缺失值，用NA指定
#------输入下面代码
cov(nudat)
colMeans(nudat)
饱和模型
首先我们构建一个饱和模型，并对MZ、DZ的每个组（twin1和twin2）估计模型的均值方差和协方差，这样我们可以得到一个基线的似然度估计，并且这个基线将会作为后续我们比较其他模型的基准。饱和模型中有1）两个均值期望的矩阵（一个是MZ的一个是DZ的），2）两个协方差期望的矩阵（即Cholesky分解），通过一个下三角矩阵乘以自己的转置。刚刚生成的双生子数据接下来通过mxData指令进行指定，模型拟合优化的幻术通过mxFIMLObjective进行指令。下面进行自模型的建模，以MZ为例。
#--------------------------------------------------------------------------------------------
# MZ submodel
#--------------------------------------------------------------------------------------------
mxModel("MZ",
#------第一个矩阵：MZ的均值期望
mxMatrix(
type="Full",
nrow=1,
ncol=2,
free=TRUE,
values=c(0,0),
name="expMeanMZ"),
#------第二个矩阵：MZ的Cholesky分解
mxMatrix(
type="Lower",
nrow=2,
ncol=2,
free=TRUE,
values=.5,
name="CholMZ"),
#------第三个矩阵：MZ的协方差期望矩阵
mxAlgebra(
CholMZ %*% t(CholMZ),
name="expCovMZ"),
mxData(
myDataMZ,
type="raw"),
mxFIMLObjective("expCovMZ", "expMeanMZ")) # 拟合函数
#--------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------
# DZ submodel, illustrating compact formatting
#--------------------------------------------------------------------------------------------
mxModel("DZ",
mxMatrix("Full", 1, 2, T, c(0,0), name="expMeanDZ"),
mxMatrix("Lower", 2, 2, T, .5, name="CholDZ"),
mxAlgebra(CholDZ %*% t(CholDZ), name="expCovDZ"),
mxData(myDataDZ, type="raw"),
mxFIMLObjective("expCovDZ", "expMeanDZ", mylabels))
#--------------------------------------------------------------------------------------------
超级模型的建立
将MZ和DZ结合起来就可以建立“super”模型，为了得到MZ和DZ数据的整体拟合似然度，我们将MZ和DZ模型的似然度的log值简单相加。超级模型里面将会包含进行这个操作的指令，这个模型叫做‘twin’。
超级的拟合
这里将对上面构建的模型进行拟合，首先是饱和模型。饱和模型中忽略双生子对之间的关系，也就是说模型中MZ和DZ双生子对之间是相同的。模型中没有潜变量，只有观测的均值和协方差。
#---------------------------------------------------------------------------------------------
# Saturated_twin_model
# DVM Bishop, 12th March 2010, based on OpenMxUsersGuide, p. 16
#---------------------------------------------------------------------------------------------
require(OpenMx)
mytwindata=read.table("mytwinfile")
#read in previously saved data created with Twin Simulate script
myDataMZ=mytwindata[,1:2] #columns 1-2 are MZ twin1 and twin2
myDataDZ=mytwindata[,3:4] #columns 3-4 are DZ twin 1 and twin 2
colnames(myDataMZ)=c("twin1","twin2")
colnames(myDataDZ)=colnames(myDataMZ)
mylabels=c("twin1","twin2")
# 模型设置Model specification starts here
mytwinSatModel <- mxModel("twinSat", # 饱和模型的配置
mxModel("MZ", #构建MZ模型
mxMatrix(type="Full",
nrow=1,
ncol= 2,
free=TRUE,
values=c(0,0),
name="expMeanMZ"), #模型中的MZ均值期望
mxMatrix("Lower",
nrow= 2,
ncol=2,
free=TRUE,
values=.5,
name="CholMZ"), #模型中的MZ的Cholesky分解
mxAlgebra(CholMZ %*% t(CholMZ), name="expCovMZ"),
mxData(myDataMZ, type="raw"),
mxFIMLObjective("expCovMZ", "expMeanMZ", mylabels)),
mxModel("DZ",
mxMatrix(type="Full",
nrow=1,
ncol= 2,
free=TRUE,
values=c(0,0),
name="expMeanDZ"),
mxMatrix(type="Lower",
nrow=2,
ncol=2,
free=TRUE,
values=.5,
name="CholDZ"),
mxAlgebra(CholDZ %*% t(CholDZ), name="expCovDZ"),
mxData(myDataDZ, type="raw"),
mxFIMLObjective("expCovDZ", "expMeanDZ", mylabels)),
#以上为最DZ建模
mxAlgebra(MZ.objective + DZ.objective, name="twin"),
# 将MZ和DZ的似然度相加
mxAlgebraObjective("twin"))
# 对mxAlgebra的表达进行估计，连个submodel同时进行估计
#---------------------------------------------------------------------------------------------------------------------------
mytwinSatFit <- mxRun(mytwinSatModel) #mxRun指令对模型进行估计
myExpMeanMZ <- mxEval(MZ.expMeanMZ, mytwinSatFit)
#估计完后，可以查看模型中参数的结果
myExpCovMZ <- mxEval(MZ.expCovMZ, mytwinSatFit)
myExpMeanDZ <- mxEval(DZ.expMeanDZ, mytwinSatFit)
myExpCovDZ <- mxEval(DZ.expCovDZ, mytwinSatFit)
LL_Sat <- mxEval(objective, mytwinSatFit)
summary(mxRun(mytwinSatModel))
#--------------------------------------------------------------------------------------------------------------------------
# compute DF for this model - this is a clunky way to do it!
msize=nrow(myDataMZ)*ncol(myDataMZ)
dsize=nrow(myDataDZ)*ncol(myDataDZ)
myDF_Sat=msize+dsize-nrow(mytwinSatFit@output$standardErrors)
#-------------------------------------------------------------------------------------------------------------------------
这里自由度可以反应样本的大小减去腰估计的变量个数。上述例子中，总的观测样本书为4000（1000对MZ和1000对DZ），要估计的参数为10个，所以自由度为3990。
- 我们没有采用RAM方法，那么需要手动画出饱和模型的路径图
- 十个要估计的参数都是啥？
- 怎么用R指令得到DZ的协方差矩阵？（只利用输出结果中DZ.CholDZ)
- Chi方没有给出数值而是NA，为什么？
第三个问题最容易回答，只需要将DZ.CholDZ中的三个值，变形成2x2的矩阵，便可以得到一个上三角：
a=matrix(c(.9968,.5562,0,.80269),nrow=2)
然后将得到的矩阵乘以自身转置，就可以得到估计的协方差矩阵。这时可以和真实值进行对比：
cod(myDataDZ)
下面的路径图便是饱和模型的路径图：
图中从“one”指向四个方框变量的路径为四个估计的均值；自环双头箭头路径指的是变量自身的方差；MZ、DZ内部的双头箭头路径指的是协方差。模型中不包含潜变量，只有观测变量（方框）和常量（三角）。**在对ACE模型进行拟合之前，我们会先检验双生子模型中的假设是否成立，即，twin1和twin2的均值和方差是否相等，并且MZ和DZ之间均值和方差是否相等。
OpenMx中一个特点是，如果两个估计的参数有相同的label，那么这两个参数将会被认为是相同的。下面对于这一点进行演示，在之前的脚本中，在expMeanMZ和expMeanDZ矩阵进行定义的位置把以下的语句加入到定义函数的name选项后：
labels="mean"
在重新跑脚本之前，先把-2LL和DF的值记录下来，并且观察重新运行脚本之后估计值的变化。我们把所有的4个均值都命名为同样的label，所以这四个均值将有相同的估计值。对于-2LL和自由度来说，会有什么改变呢？原来的脚本结果-2LL值为9911.62，DF为3990；新的脚本结果9913.154，DF为3993.通过-2LL的差和DF的差就可以决定样本的均值是否相同。
