美文网首页R语言
1.4-双生子数据中结构方程模型的应用

1.4-双生子数据中结构方程模型的应用

作者: 赵晨西西西西 | 来源:发表于2018-02-19 10:40 被阅读0次

    来源网站:http://dbtemp.blogspot.com/2011/08/sem-approach-to-twin-data.html


    对于行为遗传学分析来说,我们希望能够衡量一个比例,即人与人之间性状的相似性中有多少可以用基因的相似性来解释。在最简单的双生子分析中,我们的观测变量包括两个协方差矩阵:所有MZ同卵双胞胎一个、所有DZ异卵双胞胎一个。
    针对某一性状(变量)而言,协方差矩阵可以表示双胞胎对之间该性状的协方差(可以理解为非标准化的相关系数,反应有单位的共变程度)。
    

    这个简单模型中最重要的假设是:双胞胎之间的形状的相似性(差异)可以用遗传、共有环境(包括出生前胎内环境、家庭环境、社会经济地位、住宅区;使得双胞胎彼此更像)和他们特有的环境(这些环境因素不共有,包括一些特殊的经历[疾病、受伤等]、随机的生物效应、自身对感知的差别、以及测量的误差;造成双胞胎之间的差异)来解释。
    一个容易混淆的点是,当我们谈论个体之间遗传相似性的时候,我们认为:同卵双胞胎MZ具有完全相同遗传物质,而异卵双胞胎DZ之间有50%的相同遗传物质。这其实很难理解的(确实),因为我们也常常听到有人说我们人类和大猩猩基因相似性高达98%。这里不如这么理解:在行为遗传学分析中,我们只关注那些在不同人上有不同等位基因形式的基因,这些基因可能可以解释我看观察得到的个体差异。这些具有多态性的基因仅仅占人类整个基因组中很小的一部分。
    如果我们暂时只关心加性基因的效应(additive genetic influence),则对于异卵双胞胎可以表示如下图所示的标准ACE模型。

    标准ACE模型:DZ双胞胎
    如图所示,每个潜变量都有相对应的路径,我们的目标便是估计出这些路径。这里需要注意的是,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的差就可以决定样本的均值是否相同。

    相关文章

      网友评论

        本文标题:1.4-双生子数据中结构方程模型的应用

        本文链接:https://www.haomeiwen.com/subject/ftyqhxtx.html