精华内容
下载资源
问答
  • 随机效应与混合效应模型 SAS实践

    千次阅读 2020-08-07 09:21:13
    随机效应与混合效应模型 SAS实践两个处理变量都具有固定效应两个处理变量都具有随机效应GLM procedureMixed procedure, type 1 methodMixed procedure, REML method只有一个处理变量具有随机效应GLM ...

    精密仪器制造非常依赖零部件尺寸的精确度。一位工程师随机选取了三种仪器对某零部件的二十个关键部位进行了测量,每一个部位重复测量了两次。这个测量试验的目的是研究零部件部位与测量仪器对测量结果的影响,零部件部位与测量仪器是处理变量。实际上这两个处理变量都应该考虑为随机的,但为了比较固定效应与随机效应的区别,下面我们分别用固定效应、随机效应与混合效应模型分析这个试验。

    两个处理变量都具有固定效应

    假设两个处理变量都是固定的。

    proc glm data=randr;/*数据表读取比较长,放在文末,数据表名为randr*/
    class operator part;/*定义仪器与部位为类型变量*/
    model resp=operator|part;/*|表示使用full model,即operator part operator*part*/
    run;
    

    这个GLM模型是显著的,下面是处理变量的Type III 平方和分解:
    在这里插入图片描述
    这个结果说明测量仪器不会对测量结果造成显著的影响,但不管用什么仪器测量,每个部位测量精确度都有显著差异。

    两个处理变量都具有随机效应

    GLM procedure

    假设两个处理变量都是随机的。

    proc glm data=randr;
    class operator part;
    model resp=operator|part;
    random operator part operator*part/test; 
    run;
    

    与固定效应模型相比,随机效应模型只需要在GLM procedure中加入random语句。在这段SAS代码中,random operator part operator*part/test; 语句表示operator part operator*part这三个因素都具有随机效应,/test表示对随机效应进行F检验,之所以要加上/test是因为随机效应模型中并非所有的均方和都与残差的均方和比较:
    在这里插入图片描述
    随机效应模型也是显著的,下面是处理变量的Type III 平方和分解:
    在这里插入图片描述

    这个结果与随机模型的结果一模一样,但是我们只能用这个结果做交互项operator*part的推断,要对另外两个因素做推断,需要random语句的结果:
    在这里插入图片描述
    这两张表是random语句的结果。在第二张表中,交互项与残差的平方和分解结果与固定效应模型和随机效应模型的都相同;第一张表是random语句中/test的结果,对于operator与part而言,error term是operator*part,这个结果同样说明测量仪器不会对测量结果造成显著的影响,但每个部位测量精确度都有显著差异。

    需要注意的是random语句中的/test语句可以用test语句代替:

    proc glm data=randr;
    class operator part;
    model resp=operator|part;
    random operator part operator*part; 
    test H=operator E=operator*part;
    test H=part E=operator*part;
    run;
    

    语句test H=operator E=operator*part;表示检验operator的效应,用operator*part的均方和作为error term的均方和。语句test H=part E=operator*part;表示检验part的效应,用operator*part的均方和作为error term的均方和。这两个语句输出的结果与上面的第一张表相同:
    在这里插入图片描述

    Mixed procedure, type 1 method

    随机效应模型也可以用Mixed procedure来做,具体有两种方法,type 1 method和REML method:

    proc mixed cl maxiter=20 covtest method=type1 data = randr;/*covtest表示估计每个效应的方差、cl表示列出方差的置信区间*/
    class operator part;
    model resp = ;/*mixed procedure后的随机因素不用写出来*/
    random operator part operator*part;/*mixed procedure没有/test语句*/
    run;
    

    在这里插入图片描述
    type 1 method输出的结果与GLM procedure的random语句输出结果一样,下面是covtest与cl语句的输出,包括每个因素的方差估计、标准差以及置信区间,type 1 method估计方差的方法是最小二乘法。
    在这里插入图片描述

    Mixed procedure, REML method

    proc mixed cl maxiter=20 covtest method=reml;
    class operator part;
    model resp = ;
    random operator part operator*part;
    run;
    

    在这里插入图片描述
    REML method估计方差的方法是残差最大似然方法,是Mixed procedure的默认方法。

    只有一个处理变量具有随机效应

    假设operator是固定的因素。

    GLM procedure

    Restricted Model

    proc glm data=randr;
    class operator part;
    model resp=operator|part;
    random  part operator*part; /*operator是固定的,随机的只有part和交互项*/
    test H=operator E=operator*part;
    run;
    

    注意restricted model在做推断时,检验如下
    在这里插入图片描述
    因此上面那段SAS代码中加了一个test语句,test H=operator E=operator*part;用以对operator的效应做推断。
    在这里插入图片描述
    这是GLM procedure的输出结果,part和operator的结果是可以参考的,operator的结果需要参考test语句的输出:
    在这里插入图片描述

    Unrestricted Model

    proc glm data=randr;
    class operator part;
    model resp=operator|part;
    random  part operator*part/test; 
    run;
    

    Restricted Model与Unrestricted Model相比唯一的区别就是交互项完全随机还是交互项关于固定效应的因子的边际效应是固定的。根据Unrestricted model的EMS,
    在这里插入图片描述
    part与operator的均方和都需要与它们的交互项进行比较,所以用/test语句来实现,下面是/test的输出结果
    在这里插入图片描述

    两两比较

    proc glm;
    class operator part;
    model resp=operator|part;
    random part operator*part / test;
    means operator / tukey lines E=operator*part;
    lsmeans operator / adjust=tukey E=operator*part tdiff stderr;
    run;
    

    一般对固定效应的因素,不同的因素取值对测量值的效应两两之间可以进行比较。means语句和lsmeans语句都可以做比较,means语句means operator / tukey lines E=operator*part;表示对operator的效应进行比较,用Tukey方法,设定error term为operator*part,lines表示用线图表示结果:
    在这里插入图片描述
    Estimate是每一种测量仪器测量结果的均值,蓝线贯穿了三个均值表示这三个均值没有显著差异。lsmeans语句lsmeans operator / adjust=tukey E=operator*part tdiff stderr;可以做更精细的比较:
    在这里插入图片描述
    在这里插入图片描述

    Mixed procedure

    proc mixed alpha=.05 cl covtest;
    class operator part;
    model resp=operator / ddfm=kr;
    random part operator*part;
    lsmeans operator / alpha=.05 cl diff adjust=tukey;
    run;
    

    上面的SAS代码中只有ddfm=kr语句需要注意一下,ddfm或者ddf用来指定交互项自由度的计算方法,一般用kr就可以了。lsmeans语句的输出包括operator均值的估计、置信区间以及两两比较:
    在这里插入图片描述

    展开全文
  • 为了提高已有模型中仅单独考虑学习和(或)恶化效应所产生处理时间的准确性,提出了一个同时考虑学习和(或)恶化效应模型.该模型将学习和恶化效应函数都建模为关于位置和累积时间的函数,因此该模型较已有模型应用更...
  • 单因素试验固定效应模型方差分析 观测值的线性模型 平方和与自由度分解 例题与SPSS求解 非平衡单因素试验SPSS求解 一、观测值的线性模型 单因素试验线性可加模型为: Yij为第i个处理的第j个观测值;U为所有观测值...

    单因素试验固定效应模型方差分析

    1. 观测值的线性模型
    2. 平方和与自由度分解
    3. 例题与SPSS求解
    4. 非平衡单因素试验SPSS求解

    一、观测值的线性模型
    单因素试验线性可加模型为:

    在这里插入图片描述
    Yij为第i个处理的第j个观测值;U为所有观测值的平均值;Ti为第i个处理效应;Eij为随机误差。

    二、平方和与自由度分解
    平凡和的分解这里不做介绍,因为单因素方差分析大部分统计软件都能够计算,没有必要进行手动计算;另外还有一个原因就是CSDN编辑公式比较麻烦,这里就不敲了。PS:如果大家在CSDN有好的公式展示方法,后面会补上!~~目前需要的童鞋,自行查询相应书籍。

    自由度分解:
    自由度分解比较简单,下面详细介绍一下,大家可以根据自由度判断使用的软件给出的方差分析表是否正确,这是一个比较快速的方法。
    如果一个单因素试验,一共有a个处理,每个处理重复n次,则一共有an个观测值,具有一个约束条件所有的观测值减去观测平均值后,求和结果为0,因此总自由度为an-1;a个处理之间,每个处理观测值的平均值,和所有观测平均值之差求和为0,所以处理间自由度为a-1;而每个处理内部,n个观测值与该处理下观测平均值之差求和为0,一共有a个处理,则处理内自由度为a*(n-1)。
    总的来说,就是总自由度为观测值个数减一,也就是an-1;处理间自由度为处理个数减一,a-1;处理内自由度为a(n-1)。

    三、例题与SPSS求解
    某个树种5个不同种源种子的重量比较,分别称取100粒不同种源种子的重量,每个种源称量8次(每一批种子不同),观测值如下:
    在这里插入图片描述
    SPSS求解:
    SPSS数据表格如下:
    在这里插入图片描述
    软件操作:
    在这里插入图片描述在这里插入图片描述
    结果:
    在这里插入图片描述
    组间显著性小于0.000,说明种源种子质量之间差异极显著。

    下一章会介绍利用SPSS进行单因素非平衡试验方差分析,以及多重对比

    展开全文
  • 我们已经学习了如何处理混合效应模型。本文的重点是如何建立和可视化混合效应模型的结果。 设置 本文使用数据集,用于探索草食动物种群对珊瑚覆盖的影响。 knitr::opts_chunk$set(echo = TRUE) library...

    原文链接:http://tecdat.cn/?p=20631 

    原文出处:拓端数据部落公众号

     

    我们已经学习了如何处理混合效应模型。本文的重点是如何建立和可视化 混合效应模型的结果。

    设置

    本文使用数据集,用于探索草食动物种群对珊瑚覆盖的影响。

    knitr::opts_chunk$set(echo = TRUE)
    
    library(tidyverse) # 数据处理
    library(lme4) #  lmer   glmer 模型
    
    
    
    me_data <- read_csv("mixede.csv")

    创建一个基本的混合效应模型:

    该模型以珊瑚覆盖层为因变量(elkhorn_LAI),草食动物种群和深度为固定效应(c。 urchinden,c.fishmass,c.maxD)和调查地点作为随机效应(地点)。

    注意:由于食草动物种群的测量规模存在差异,因此我们使用标准化的值,否则模型将无法收敛。我们还使用了因变量的对数。我正在根据这项特定研究对数据进行分组。

    summary(mod)
    ## Linear mixed model fit by maximum likelihood  ['lmerMod']
    
    ## 
    ##      AIC      BIC   logLik deviance df.resid 
    ##    116.3    125.1    -52.1    104.3       26 
    ## 
    ## Scaled residuals: 
    ##     Min      1Q  Median      3Q     Max 
    ## -1.7501 -0.6725 -0.1219  0.6223  1.7882 
    ## 
    ## Random effects:
    ##  Groups   Name        Variance Std.Dev.
    ##  site     (Intercept) 0.000    0.000   
    ##  Residual             1.522    1.234   
    ## Number of obs: 32, groups:  site, 9
    ## 
    ## Fixed effects:
    ##             Estimate Std. Error t value
    ## (Intercept)  10.1272     0.2670  37.929
    ## c.urchinden   0.5414     0.2303   2.351
    ## c.fishmass    0.4624     0.4090   1.130
    ## c.maxD        0.3989     0.4286   0.931
    ## 
    ## Correlation of Fixed Effects:
    ##             (Intr) c.rchn c.fshm
    ## c.urchinden  0.036              
    ## c.fishmass  -0.193  0.020       
    ## c.maxD       0.511  0.491 -0.431
    ## convergence code: 0
    ## boundary (singular) fit: see ?isSingular

    绘制效应大小图:

    如果您有很多固定效应,这很有用。

     

    plot(mod)

    效应大小的格式化图:

    让我们更改轴标签和标题。

    # 注意:轴标签应按从下到上的顺序排列。 
    # 要查看效应大小和p值,设置show.values和show.p= TRUE。只有当效应大小的值过大时,才会显示P值。
      title="草食动物对珊瑚覆盖的影响")

    模型结果表输出:

    创建模型摘要输出表。这将提供预测变量,包括其估计值,置信区间,估计值的p值以及随机效应信息。

     

     tab(mod)

    格式化表格

    # 注:预测标签(pred.labs)应从上到下排列;dv.labs位于表格顶部的因变量的名称。
     
     
                      pred.labels =c("(Intercept)", "Urchins", "Fish", "Depth"),
     

     

     

    用数据绘制模型估计

    我们可以在实际数据上绘制模型估计值!我们一次只针对一个变量执行此操作。注意:数据已标准化以便在模型中使用,因此我们绘制的是标准化数据值,而不是原始数据

    步骤1:将效应大小估算值保存到data.frame中

    # 使用函数。 term=固定效应,mod=你的模型。
    
    effect(term= "c.urchinden", mod= mod)
    summary(effects) #值的输出
    ## 
    ##  c.urchinden effect
    ## c.urchinden
    ##     -0.7      0.4        2        3        4 
    ##  9.53159 10.12715 10.99342 11.53484 12.07626 
    ## 
    ##  Lower 95 Percent Confidence Limits
    ## c.urchinden
    ##      -0.7       0.4         2         3         4 
    ##  8.857169  9.680160 10.104459 10.216537 10.306881 
    ## 
    ##  Upper 95 Percent Confidence Limits
    ## c.urchinden
    ##     -0.7      0.4        2        3        4 
    ## 10.20601 10.57414 11.88238 12.85314 13.84563
    # 将效应值另存为df:
    x  <- as.data.frame(effects)

    步骤2:使用效应值df绘制估算值

    如果要保存基本图(仅固定效应和因变量数据),可以将其分解为单独的步骤。注意:对于该图,我正在基于此特定研究对数据进行分组。

    #基本步骤:
      #1创建空图
    
      #2 从数据中添加geom_points()
    
      #3 为模型估计添加geom_point。我们改变颜色,使它们与数据区分开来
    
      #4 为MODEL的估计值添加geom_line。改变颜色以配合估计点。
    
      #5 添加具有模型估计置信区间的geom_ribbon
    
      #6 根据需要编辑标签!
    
    #1
    chin_plot <- ggplot() + 
      #2
    geom_point(data ,  + 
      #3
      geom_point(data=x_, aes(x= chinde, y=fit), color="blue") +
      #4
      geom_line(data=x, aes(x= chinde, y=fit), color="blue") +
      #5
      geom_ribbon(data= x , aes(x=c.urchinden, ymin=lower, ymax=upper), alpha= 0.3, fill="blue") +
      #6
      labs(x="海胆(标准化)", y="珊瑚覆盖层")
    
    chin_plot


    最受欢迎的见解

    1.基于R语言的lmer混合线性回归模型

    2.R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)

    3.R语言线性混合效应模型实战案例

    4.R语言线性混合效应模型实战案例2

    5.R语言线性混合效应模型实战案例

    6.线性混合效应模型Linear Mixed-Effects Models的部分折叠Gibbs采样

    7.R语言LME4混合效应模型研究教师的受欢迎程度

    8.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

    9.使用SAS,Stata,HLM,R,SPSS和Mplus的分层线性模型HLM

    展开全文
  • R语言 线性混合效应模型实战案例

    千次阅读 2019-06-13 08:47:32
    处理分组数据和复杂层次结构的分析师,从嵌入在参与者中的测量,嵌套在州内的县或嵌套在教室内的学生,经常发现他们需要建模工具来反映他们数据的这种结构。在R中,有两种主要的方法来拟合多级模型,这些模型考虑了...

    原文链接:http://tecdat.cn/?p=3059

    原文出处:拓端数据部落公众号

     

    介绍

    处理分组数据和复杂层次结构的分析师,从嵌入在参与者中的测量,嵌套在州内的县或嵌套在教室内的学生,经常发现他们需要建模工具来反映他们数据的这种结构。在R中,有两种主要的方法来拟合多级模型,这些模型考虑了数据中的这种结构。这些教程将向用户展示如何使用lme4R中的包来拟合线性和非线性混合效果模型,以及如何使用rstan以完全适合贝叶斯多级模型。这里的重点是如何使模型适合R而不是模型背后的理论。有关多级建模的背景知识,请参阅参考资料。 

    本教程将介绍如何lme4  设置和运行一些基本模型,其中包括:

    • 在R中构造变化的截距,变化的斜率以及变化的斜率和截距模型
    • 从混合效应模型中生成预测和解释参数
    • 广义和非线性多层次模型
    • 完全贝叶斯多级模型适合rstan或其他MCMC方法

    设置 环境

    在R中开始多级建模很简单。lme4是在R中实现多级模型的规范包,尽管有许多包依赖并增强其功能集,包括贝叶斯扩展。lme4 最近已被重写以提高速度并整合C ++代码库,因此封装的功能有些不断变化。 

    要安装lme4,我们只需运行:

    # 主要版本
    install.packages("lme4")
    
    # 或安装开发版本
    library(devtools)
    install_github("lme4", user = "lme4")
    

    读入数据

    多级模型适用于特定类型的数据结构,其中单元嵌套在组内(通常为5个以上组),并且我们希望对数据的组结构进行建模。对于我们的介绍性示例,我们将从lme4文档中的一个简单示例开始,并解释模型正在执行的操作。 

    library(lme4)  # 加载库
    library(arm)  # R中用于回归的函数
      # summary(lmm.data)
    head(lmm.data)
    
    ##   id extro  open agree social class school
    ## 1  1 63.69 43.43 38.03  75.06     d     IV
    ## 2  2 69.48 46.87 31.49  98.13     a     VI
    ## 3  3 79.74 32.27 40.21 116.34     d     VI
    ## 4  4 62.97 44.41 30.51  90.47     c     IV
    ## 5  5 64.25 36.86 37.44  98.52     d     IV
    ## 6  6 50.97 46.26 38.83  75.22     d      I
    

     模型

    让我们首先拟合一个简单的OLS回归 。 

    OLSexamp <- lm(extro ~ open + agree + social, data = lmm.data)
    display(OLSexamp)
    
    ## lm(formula = extro ~ open + agree + social, data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 57.84     3.15  
    ## open         0.02     0.05  
    ## agree        0.03     0.05  
    ## social       0.01     0.02  
    ## ---
    ## n = 1200, k = 4
    ## residual sd = 9.34, R-Squared = 0.00
    

     R模型接口非常简单,首先指定因变量,然后是 ~符号,预测变量,每个都被命名。加法符号表明这些被建模为加性效应。最后,我们指定要计算模型的数据。这里我们使用该lm函数执行OLS回归,但R中还有许多其他选项。

    如果我们想要提取诸如AIC之类的度量 。

    MLexamp <- glm(extro ~ open + agree + social, data = lmm.data)
    display(MLexamp)
    
    ## glm(formula = extro ~ open + agree + social, data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 57.84     3.15  
    ## open         0.02     0.05  
    ## agree        0.03     0.05  
    ## social       0.01     0.02  
    ## ---
    ##   n = 1200, k = 4
    ##   residual deviance = 104378.2, null deviance = 104432.7 (difference = 54.5)
    ##   overdispersion parameter = 87.3
    ##   residual sd is sqrt(overdispersion) = 9.34
    
    AIC(MLexamp)
    
    ## [1] 8774
    

    这导致模型拟合较差。现在让我们看一个简单的模型。

    拟合不同的 模型

    我们的下一步可能是使用分组变量(如学校或班级)来拟合不同的 模型。 

    MLexamp.2 <- glm(extro ~ open + agree + social + class, data = lmm.data)
    display(MLexamp.2)
    
    ## glm(formula = extro ~ open + agree + social + class, data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 56.05     3.09  
    ## open         0.03     0.05  
    ## agree       -0.01     0.05  
    ## social       0.01     0.02  
    ## classb       2.06     0.75  
    ## classc       3.70     0.75  
    ## classd       5.67     0.75  
    ## ---
    ##   n = 1200, k = 7
    ##   residual deviance = 99187.7, null deviance = 104432.7 (difference = 5245.0)
    ##   overdispersion parameter = 83.1
    ##   residual sd is sqrt(overdispersion) = 9.12
    
    AIC(MLexamp.2)
    
    ## [1] 8719
    
    anova(MLexamp, MLexamp.2, test = "F")
    
    ## Analysis of Deviance Table
    ## 
    ## Model 1: extro ~ open + agree + social
    ## Model 2: extro ~ open + agree + social + class
    ##   Resid. Df Resid. Dev Df Deviance    F  Pr(>F)    
    ## 1      1196     104378                             
    ## 2      1193      99188  3     5190 20.8 3.8e-13 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    

    这通常被称为固定效应 。 

    MLexamp.3 <- glm(extro ~ open + agree + social + school, data = lmm.data)
    display(MLexamp.3)
    
    ## glm(formula = extro ~ open + agree + social + school, data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 45.02     0.92  
    ## open         0.01     0.01  
    ## agree        0.03     0.02  
    ## social       0.00     0.00  
    ## schoolII     7.91     0.27  
    ## schoolIII   12.12     0.27  
    ## schoolIV    16.06     0.27  
    ## schoolV     20.43     0.27  
    ## schoolVI    28.05     0.27  
    ## ---
    ##   n = 1200, k = 9
    ##   residual deviance = 8496.2, null deviance = 104432.7 (difference = 95936.5)
    ##   overdispersion parameter = 7.1
    ##   residual sd is sqrt(overdispersion) = 2.67
    
    AIC(MLexamp.3)
    
    ## [1] 5774
    
    anova(MLexamp, MLexamp.3, test = "F")
    
    ## Analysis of Deviance Table
    ## 
    ## Model 1: extro ~ open + agree + social
    ## Model 2: extro ~ open + agree + social + school
    ##   Resid. Df Resid. Dev Df Deviance    F Pr(>F)    
    ## 1      1196     104378                            
    ## 2      1191       8496  5    95882 2688 <2e-16 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    

    学校因子大大提高了我们的模型拟合度。但是,我们如何解释这些影响呢?

    table(lmm.data$school, lmm.data$class)
    
    ##      
    ##        a  b  c  d
    ##   I   50 50 50 50
    ##   II  50 50 50 50
    ##   III 50 50 50 50
    ##   IV  50 50 50 50
    ##   V   50 50 50 50
    ##   VI  50 50 50 50
    

    在这里,我们可以看到我们有一个完美平衡的设计,在课堂和学校的每个组合中有50个观察结果 。

    让我们尝试对这些独特的因素进行建模。 

    MLexamp.4 <- glm(extro ~ open + agree + social + school:class, data = lmm.data)
    display(MLexamp.4)
    
    ## glm(formula = extro ~ open + agree + social + school:class, data = lmm.data)
    ##                  coef.est coef.se
    ## (Intercept)       80.36     0.37 
    ## open               0.01     0.00 
    ## agree             -0.01     0.01 
    ## social             0.00     0.00 
    ## schoolI:classa   -40.39     0.20 
    ## schoolII:classa  -28.15     0.20 
    ## schoolIII:classa -23.58     0.20 
    ## schoolIV:classa  -19.76     0.20 
    ## schoolV:classa   -15.50     0.20 
    ## schoolVI:classa  -10.46     0.20 
    ## schoolI:classb   -34.60     0.20 
    ## schoolII:classb  -26.76     0.20 
    ## schoolIII:classb -22.59     0.20 
    ## schoolIV:classb  -18.71     0.20 
    ## schoolV:classb   -14.31     0.20 
    ## schoolVI:classb   -8.54     0.20 
    ## schoolI:classc   -31.86     0.20 
    ## schoolII:classc  -25.64     0.20 
    ## schoolIII:classc -21.58     0.20 
    ## schoolIV:classc  -17.58     0.20 
    ## schoolV:classc   -13.38     0.20 
    ## schoolVI:classc   -5.58     0.20 
    ## schoolI:classd   -30.00     0.20 
    ## schoolII:classd  -24.57     0.20 
    ## schoolIII:classd -20.64     0.20 
    ## schoolIV:classd  -16.60     0.20 
    ## schoolV:classd   -12.04     0.20 
    ## ---
    ##   n = 1200, k = 27
    ##   residual deviance = 1135.9, null deviance = 104432.7 (difference = 103296.8)
    ##   overdispersion parameter = 1.0
    ##   residual sd is sqrt(overdispersion) = 0.98
    
    AIC(MLexamp.4)
    
    ## [1] 3396
    

    这非常有用,但如果我们想了解学校的影响和课堂的影响,以及学校和班级的影响,该怎么办? 

    MLexamp.5 <- glm(extro ~ open + agree + social + school * class - 1, data = lmm.data)
    display(MLexamp.5)
    
    ## glm(formula = extro ~ open + agree + social + school * class - 
    ##     1, data = lmm.data)
    ##                  coef.est coef.se
    ## open              0.01     0.00  
    ## agree            -0.01     0.01  
    ## social            0.00     0.00  
    ## schoolI          39.96     0.36  
    ## schoolII         52.21     0.36  
    ## schoolIII        56.78     0.36  
    ## schoolIV         60.60     0.36  
    ## schoolV          64.86     0.36  
    ## schoolVI         69.90     0.36  
    ## classb            5.79     0.20  
    ## classc            8.53     0.20  
    ## classd           10.39     0.20  
    ## schoolII:classb  -4.40     0.28  
    ## schoolIII:classb -4.80     0.28  
    ## schoolIV:classb  -4.74     0.28  
    ## schoolV:classb   -4.60     0.28  
    ## schoolVI:classb  -3.87     0.28  
    ## schoolII:classc  -6.02     0.28  
    ## schoolIII:classc -6.54     0.28  
    ## schoolIV:classc  -6.36     0.28  
    ## schoolV:classc   -6.41     0.28  
    ## schoolVI:classc  -3.65     0.28  
    ## schoolII:classd  -6.81     0.28  
    ## schoolIII:classd -7.45     0.28  
    ## schoolIV:classd  -7.24     0.28  
    ## schoolV:classd   -6.93     0.28  
    ## schoolVI:classd   0.06     0.28  
    ## ---
    ##   n = 1200, k = 27
    ##   residual deviance = 1135.9, null deviance = 4463029.9 (difference = 4461894.0)
    ##   overdispersion parameter = 1.0
    ##   residual sd is sqrt(overdispersion) = 0.98
    
    AIC(MLexamp.5)
    
    ## [1] 3396
    

    探索随机斜率

    另一种选择是为每个学校和班级组合建立一个单独的模型。如果我们认为我们的变量之间的关系可能高度依赖于学校和班级组合,我们可以简单地拟合一系列模型并探索它们之间的参数变化:

    require(plyr)
    
      display(modellist[[1]])
    
    ## glm(formula = extro ~ open + agree + social, data = x)
    ##             coef.est coef.se
    ## (Intercept) 35.87     5.90  
    ## open         0.05     0.09  
    ## agree        0.02     0.10  
    ## social       0.01     0.03  
    ## ---
    ##   n = 50, k = 4
    ##   residual deviance = 500.2, null deviance = 506.2 (difference = 5.9)
    ##   overdispersion parameter = 10.9
    ##   residual sd is sqrt(overdispersion) = 3.30
    
    display(modellist[[2]])
    
    ## glm(formula = extro ~ open + agree + social, data = x)
    ##             coef.est coef.se
    ## (Intercept) 47.96     2.16  
    ## open        -0.01     0.03  
    ## agree       -0.03     0.03  
    ## social      -0.01     0.01  
    ## ---
    ##   n = 50, k = 4
    ##   residual deviance = 47.9, null deviance = 49.1 (difference = 1.2)
    ##   overdispersion parameter = 1.0
    ##   residual sd is sqrt(overdispersion) = 1.02
    

    我们将在未来的教程中更深入地讨论此策略,包括如何在此命令中生成的模型列表中进行性能推断。

    建立不同的斜率模型

    虽然上述所有技术都是解决这一问题的有效方法,但当我们明确感兴趣的是群体之间的变化时,它们并不一定是最好的方法。这是混合效果建模框架有用的地方。现在我们使用lmer具有熟悉的公式接口的函数, 使用特殊语法指定组级变量:(1|school) ,使lmer拟合具有变量截距组效果的线性模型school

     display(MLexamp.6)
    
    ## lmer(formula = extro ~ open + agree + social + (1 | school), 
    ##     data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 59.12     4.10  
    ## open         0.01     0.01  
    ## agree        0.03     0.02  
    ## social       0.00     0.00  
    ## 
    ## Error terms:
    ##  Groups   Name        Std.Dev.
    ##  school   (Intercept) 9.79    
    ##  Residual             2.67    
    ## ---
    ## number of obs: 1200, groups: school, 6
    ## AIC = 5836.1, DIC = 5789
    ## deviance = 5806.5
    

    我们可以使用多个组来拟合多个组效果。

     display(MLexamp.7)
    
    ## lmer(formula = extro ~ open + agree + social + (1 | school) + 
    ##     (1 | class), data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 60.20     4.21  
    ## open         0.01     0.01  
    ## agree       -0.01     0.01  
    ## social       0.00     0.00  
    ## 
    ## Error terms:
    ##  Groups   Name        Std.Dev.
    ##  school   (Intercept) 9.79    
    ##  class    (Intercept) 2.41    
    ##  Residual             1.67    
    ## ---
    ## number of obs: 1200, groups: school, 6; class, 4
    ## AIC = 4737.9, DIC = 4683
    ## deviance = 4703.6
    

    最后,我们可以通过以下语法拟合嵌套:

     display(MLexamp.8)
    
    ## lmer(formula = extro ~ open + agree + social + (1 | school/class), 
    ##     data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 60.24     4.01  
    ## open         0.01     0.00  
    ## agree       -0.01     0.01  
    ## social       0.00     0.00  
    ## 
    ## Error terms:
    ##  Groups       Name        Std.Dev.
    ##  class:school (Intercept) 2.86    
    ##  school       (Intercept) 9.69    
    ##  Residual                 0.98    
    ## ---
    ## number of obs: 1200, groups: class:school, 24; school, 6
    ## AIC = 3568.6, DIC = 3508
    ## deviance = 3531.1
    
    
    

    在这里(1|school/class),我们想要为1|学校和学校内的课程设置不同截距的混合效应。

    用lmer拟合变化的斜率模型

    但是,如果我们想要探索不同学生水平指标的影响,因为它们因教室而异。我们可以拟合不同的斜率模型,而不是按学校(或学校/班级)拟合模型。在这里,我们修改我们的随机效应项,在分组术语之前包含变量:(1 + open|school/class)告诉R拟合变化的斜率和不同的学校和学校类别的截距模型,并允许open变量的斜率因学校而异。

     
    display(MLexamp.9)
    
    ## lmer(formula = extro ~ open + agree + social + (1 + open | school/class), 
    ##     data = lmm.data)
    ##             coef.est coef.se
    ## (Intercept) 60.26     3.93  
    ## open         0.01     0.01  
    ## agree       -0.01     0.01  
    ## social       0.00     0.00  
    ## 
    ## Error terms:
    ##  Groups       Name        Std.Dev. Corr 
    ##  class:school (Intercept) 2.62          
    ##               open        0.01     1.00 
    ##  school       (Intercept) 9.51          
    ##               open        0.00     1.00 
    ##  Residual                 0.98          
    ## ---
    ## number of obs: 1200, groups: class:school, 24; school, 6
    ## AIC = 3574.7, DIC = 3506
    ## deviance = 3529.3
    

    结论

    在R语言和生态系统中,拟合混合效应模型和探索组变异非常容易。在以后的教程中,我们将探索模型的比较,使用混合效果模型进行推理,以及创建混合效果模型的图形表示了解它们的效果。

    附录

    
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## 
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base     
    ## 
    ## other attached packages:
    ## [1] plyr_1.8        arm_1.6-10      MASS_7.3-29     lme4_1.0-5     
    ## [5] Matrix_1.1-0    lattice_0.20-24 knitr_1.5      
    ## 
    ## loaded via a namespace (and not attached):
    ##  [1] abind_1.4-0    coda_0.16-1    evaluate_0.5.1 formatR_0.10  
    ##  [5] grid_3.0.1     minqa_1.2.1    nlme_3.1-113   splines_3.0.1 
    ##  [9] stringr_0.6.2  tools_3.0.1

     

     

    非常感谢您阅读本文,有任何问题请在下面留言!


    参考文献

     

    1.基于R语言的lmer混合线性回归模型

     

    2.R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)

     

    3.R语言线性混合效应模型实战案例

     

    4.R语言线性混合效应模型实战案例2

     

    5.R语言线性混合效应模型实战案例

     

    6.线性混合效应模型Linear Mixed-Effects Models的部分折叠Gibbs采样

     

    7.R语言LME4混合效应模型研究教师的受欢迎程度

     

    8.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

     

    9.使用SAS,Stata,HLM,R,SPSS和Mplus的分层线性模型HLM

    展开全文
  • 本文介绍了线性混合效应模型的新型贝叶斯分析。该分析基于部分折叠的方法,该方法允许某些组件从模型中部分折叠。得到的部分折叠的Gibbs(PCG)采样器被构造成适合线性...已经开发出混合效应模型处理相关响应数...
  • 在本次网络研讨会中,您将学习如何使用 Simscape 在信号处理和通信系统模型的环境中模拟模拟基带电路的效果。 Simscape 使用用于建模和仿真多域物理系统(包括模拟电路的效果)的工具扩展了 Simulink:registered:。...
  • 提出了一种比较准确的仿真模型用来估计把非线性噪声当作加性高斯处理的高斯噪声(GN)模型的适用范围,同时对该模型进行了大量的仿真验证,包括不同的光纤类型和不同的调制格式。仿真结果证明,GN模型的非线性噪声功率...
  • 配对设计:同一对子的两个实验单位可以随机分配处理,两个实验单位同期观察试验结果,可以比较处理组间差别;要求同一对子的两个实验单位的观察结果分别与差值相互独立,差值服从正态分布;用平均差值推论处理的作用...
  • 在面板数据下基于变系数固定效应模型,采用两步估计方法和截面核估计方法来估计模型中的变系数部分及...本文通过两组不同样本的模拟数据验证了两种估计方法的有效性,说明变系数固定效应模型处理面板数据的有效工具。
  • 值得注意的是,这是第一个可以使用受限估计最大似然 (REML) 处理具有固定和随机效应的方差分析和方差分量模型的 Matlab 统计软件包,包括小样本量的现代校正。 这允许改进对固定(最佳线性无偏估计)和随机项(最佳...
  • 处理了包含克尔效应和斯塔克效应的原子与场相互作用模型。由于在其有效哈密顿中作了旋转波近似,两种效应的影响可以归结为对原子与场失谐的改变,这种失谐变化依赖于一个系统的守恒量。研究表明,克尔效应和斯塔克效应...
  • 来自 McKay, JA (1998) 的数值模型。 直接探测多普勒测风激光雷达的建模。 I.边缘技术。 应用光学,37(27),6480-6486。... 该模型包括缺陷和系统光扩展的处理,假设来自光源的通光Kong径的照明均匀。
  • 具体来说,我们主张将闪烁体探测器中质子上的弹性散射引起的事件与振荡效应不相关,该事件对振荡效应不敏感,可以用作模型独立的归一化处理,该事件应与νe的反beta衰减引起的事件进行比较。 在水切伦科夫探测器...
  • 讨论了参数效应非线性贡献因子的处理方法.基于真误差的定义,给出了GMSE的定义,以此为基础,提出了以新的平差准则--GSMSE为基础的非线性平差模型的综合处理方法.作为实例,在非线性模型空间,借助迭代技术,研究了...
  • 在方差随机效应模型的单向分析中,人们会期望测量值靠得很近,因为测量的是相同的数量,并指出类效应源于共同的来源,因此不应有太大差异。 在随机模型中,允许存在一些变化,并且通过假设随机类效应的分布来辅助...
  • 本文介绍了场效应晶体管在低电平下变阻特性的物理概念及其用途;并叙述了如何通过实验,测出实验数据,然后根据数据处理原理并借助于计算机求出其数学模型
  • 在这样的非线性模型中,未观察变量的自相关导致包含高维积分的难以处理的似然。 为了解决这个问题,我们使用涉及低得多的积分阶的复合似然。 然而,参数识别变得有问题,因为在低维分布中使用的信息可能不够丰富以...
  • 在GM(1,1)模型预测中,缓冲算子能有效处理含有冲击扰动因素的原始数据,改善模型的预测效果. 本文在系统分析缓冲算子对GM(1,1)预测作用过程的基础上,提出了GM(1...
  • 曲面被离散成许多的平面多边形,如果我们的网格较大,离散度较粗,在模型表面使用明暗处理后,两两相邻的多边形会出现凸起或者是凹陷的折痕,在连接处显得比周围处亮或者暗,这就是所谓的马赫夫效应,如下图所示 ...
  • 建立强可忽略处理分配条件下因果推断的结构回归模型,估计平均处理效应.在正态分布假设下,总体参数的极大似然估计是渐近正态无偏估计.提出的方法可推广到具有测量误差的一般情形.
  • 道愿望产生的时间累积效应模型,克服了传统换道模型只根据瞬时速度利益进 行换道决策的不合理性。 4. 基于神经网络的跟驰换道模型研究。跟驰过程中驾驶员所感知的信息与 驾驶员反应之间的映射关系是非线性的,采用...
  • 文中在基于VAR模型下,将技术进步进一步分解为人力资本、R&D支出以及技术市场成交额3个因素,探讨技术进步与我国经济增长之间的动态相关关系,有别于一般选取技术进步的单一因素作为研究对象来说明与经济增长效应之间的...
  • 提出了一种基于增广参数Kalman滤波的多路径效应系统误差估计方法,将系统误差作为状态参数,并对其建立一阶AR模型,同时利用多路径重复性特性,更新多路径误差改正模型,在一定程度上解决了固定多路径误差模型随着...
  • 利用 MATLAB 仿真多普勒效应 某某某 摘 要分析多普勒效应特性建立数学模型利用 MATLAB 软件对其进行仿真试验进 行定量分析根据仿真试验结果绘制出听者接收到的信号的频率变化曲线以及用信号处理 工具箱函数 ...
  • 平稳性检验方法:时序图检验、自相关图、偏自相关图检验、adfuller...非平稳序列确定性分析:趋势分析、季节效应(周期性)分析、综合分析、X11分析 非平稳序列转换为平稳序列:差分、方差齐性变换、平滑、分解 ...
  • 振铃效应

    2016-10-26 17:42:38
    是由于在图像复原中选取了不适当的图像模型造成的,振铃效应产生的直接原因是图像退化过程中信息量的丢失,尤其是高频信息的丢失,其严重降低了复原图像的质量,并且使得难于对复原图像进行后续处理。 振铃...
  • 版权声明:本文为博主原创文章,未经博主允许不得转载。...2.4 语音的感知  2.4.1 几个概念  语音的听觉感知是一个复杂的人脑...听觉感知的试验主要还在测试响度、音高和掩蔽效应等。人耳听觉界限的范围大约为20Hz~...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 545
精华内容 218
关键字:

处理效应模型