加权变换可以消除异方差性，使随机误差项变成同方差的。这样才会满足线性回归模型的经典假定——同方差性。因为加权最小二乘将待估计参数的估计方差的各分量变为了不相关的。
来源
模型中缺少某些解释变量，从而随机扰动项产生系统模式，由于随机扰动项ui包含了所有无法用解释变量表示的各种因素对被解释变量的影响，即模型中略去的经济变量对被解释变量的影响。如果其中被略去的某一因素或某些因素随着解释变量观测值的不同而对被解释变量产生不同的影响，就会使ui产生异方差性。
具体回答如图：
当一个随机二维向量的两个分量呈独立的、有着相同的方差的正态分布时，这个向量的模呈瑞利分布。
扩展资料：
如果变量可以在某个区间内取任一实数，即变量的取值可以是连续的，这随机变量就称为连续型随机变量。例如，公共汽车每15分钟一班，某人在站台等车时间x是个随机变量，x的取值范围是[0,15），它是一个区间，从理论上说在这个区间内可取任一实数35等。
对于连续型随机变量X，若其定义域为（a，b），概率密度函数为f（x），连续型随机变量X方差计算公式：D（X）=（x-μ）^2 f（x） dx
方差刻画了随机变量的取值对于其数学期望的离散程度。（标准差、方差越大，离散程度越大）
若X的取值比较集中，则方差D（X）较小，若X的取值比较分散，则方差D（X）较大。
因此，D（X）是刻画X取值分散程度的一个量，它是衡量取值分散程度的一个尺度。
参考资料来源：百度百科--瑞利分布
L(θ)=f(X1)f(X2)f(X3)f(Xn)=6^n/θ^(2n)Π(θ-xi)
l(θ)=ln[L(θ)]=nln6-2nlnθ+Σln(θ-xi)
0=dl(θ)/dθ=-2n/θ+Σ1/(θ-xi)
主成分分析的主要思想是将样本数据投影到一个维bai数较低的正交子空间内，而投影后的数据又能尽可能多的表达原来数据的波动情况（方差）
对于一个线性变换duA，成立Var(Ax)=AVar(x)A^T
设变量x的协方差矩阵为M。M为对称半正定矩阵，可以对角化 M=QDQ^dao-1，其中Q是正交矩阵，D是对焦矩阵。
如果选取正交变换y=Q^-1x，根据上面给出的方差公式，变换后版的数据方差为Q^-1QDQ^-1Q=D是一个对角矩阵（设D的元素从上到下递减排列），其方差为对角线上的元素即原变量协方差矩阵M的特征值，实际做的时候要舍去方差较小的几个维度。
扩展资料：
主成分分析作为基础的数学分析方法，其实际应用十分广泛，比如人口统计学、数量地理学、分子动力学模拟、数学建模、数理分析等学科中均有应用，是一种常用的多变量分析方法。
主成分分析是设法将原来众多具有一定相关性（比如P个指标），重新组合成一组新的互相无关的综合指标来代替原来的指标。
主成分分析，是考察多个变量间相关性一种多元统计方法，研究如何通过少数几个主成分来揭示多个变量间的内部结构，即从原始变量中导出少数几个主成分，使它们尽可能多地保留原始变量的信息，且彼此间互不相关通常数学上的处理就是将原来P个指标作线性组合，作为新的综合指标。
参考资料来源：百度百科-主成分分析
有线性混合模型简介
混合线性模型(mixed linear model)是一种方差分量模型。在方差分量模型中，把既含有固定效应，又含有随机效应的模型，称为混合线性模型 信息来源：百度 一般线性模型中仅包含固定效应和噪声两项影响因素，也就是不会考虑随机因素对模型的影响。 数据集：R的MASS包中oast数据集 数据解释 ： 通过三个品种四个水平的施肥处理，对燕麦进行了小区试验。实验分为6个区块，每个区块分为4个子区块。主小区采用品种，次小区采用施肥处理。 线性模型是需要求出关于在因变量上的权重(斜率)，那么针对分类型的因变量要怎么计算呢？——创建虚拟变量代替我们的分类水平进行运算。
虚拟变量
因变量属于分类变量，所以在建立模型之前需要创建虚拟变量：
虚拟变量的作用就是对固定因子的因子进行量化。例如：在varieties中有三种不同的品种分别为：Goldenrain、Marvellous、Victory，可以通过R的contrasts()函数查看或设置不同品种的虚拟变量。
默认采用treatment编码方式，每一列的求和等于1，至于为什么Goldenrain该行全都是0，因为默认情况下，第一行是其他组的参照组。注：若是采用sum编码方式：每一列的求和等于0。
展示一下treatment编码方式和sum编码方式。注：contrsum是定义无序因子的，若为有序可以采用contrpoly来定义。
例如：在3水平的因子变量上展示
同理nitrogen变量的虚拟变量编码方式如下：
观察两张可以得到关于虚拟变量创建的规律：
1、对于某一个因子型变量包含K个水平的因子，那么会产生K-1个对比向量；
2、对比向量以一个K行(K-1)列的矩阵组成。
3、默认的情况下以第一行作为其他因子的参照。(比如：Goldenrain，00cwt)
那么有了虚拟变量后要如何计算呢？
1、单因素多水平(treatment编码方式)：
运算公式：y=b+wx+beta ;其中，b是截距，w是权重(斜率),beta(随机误差)
Y = b+w1x1+w2x1+beta
A水平：YA = b+w10+w2+0+beta = b+beta
B水平：YB = b+w11+w20+beta = b+w1+beta
C水平：YC = b+w10+w21+beta = b+w2+beta
那么模型的回归系数就等于：
W1=YB — YA
W2=YC — YA
结论：回归系数等于真实效应差异。
2、两因素两水平有交互作用的情况(sum编码方式)：
Y=b+w1Xa+w2Xb+w3XaXb+beta
主效应的回归系数：
Ya1=b-w1-w2+beta
Ya2=b+w1-w2+beta
W1=(Ya2-Ya1)/2
同理：
W2=(Yb2-Yb1)/2
结论：主效应的回归系数等于真实效应的一半。
交互效应的回归系数：
Ya1b1=b-w1-w2+w3+beta
Ya1b2=b-w1+w2-w3+beta
Ya2b1=b+w1-w2-w3+beta
Ya2b2=b+w1+w2+w3+beta
W3=((Ya2b2-Ya2b1)-(Ya1b2-Ya1b1))/4
结论：交互效应等于真实交互效应的四分之一
你们也可以自己去尝试使用treatment编码方式，也就是将上面的Xa、Xb中的-1改为0。若不考虑交互作用的话，将公式中的w3XaXb剔除，再运算即可。
treatment编码方式做两因素两水平无交互作用时，其主效应的回归系数是等于真实效应的。若有交互作用时，其主效应的回归系数是不等于真实效应的，但交互作用的回归系数等于真实效应。
那么对于本案例中属于双因素多水平的怎么计算回归系数呢？
3、K因素N水平的情况下，如何保证回归系数等于真实的效应值？
R里面是没有这个函数的，可以自己设置编码。即每个因素产生N(N-1)的矩阵，其中每一列代表一个比较，每一行代表一个水平。第一个水平一般是作为为参照组，其编码为-(1/N)，每一列对应基线之外的某个水平与基线的对比，所在的行编码为1-(1/N)，其他行编码为-(1/N)。你们可以带入公式去验证一下。
DV m for(i in 1:n){ for(j in 1:n-1){ if(i-1 == j) m[i,j] else m[i,j] } } return(m)} DV(3) [,1] [,2][1,] -03333333 -03333333[2,] 06666667 -03333333[3,] -03333333 06666667DV(4) [,1] [,2] [,3][1,] -025 -025 -025[2,] 075 -025 -025[3,] -025 075 -025[4,] -025 -025 075#将虚拟变量赋值给案例中固定因子contrasts(oats$varieties) contrasts(oats$nitrogen) contrasts(oats$varieties) [,1] [,2]Goldenrain -03333333 -03333333Marvellous 06666667 -03333333Victory -03333333 06666667contrasts(oats$nitrogen) [,1] [,2] [,3]00cwt -025 -025 -02502cwt 075 -025 -02504cwt -025 075 -02506cwt -025 -025 075
其中我们将blocks变量作为被试，那么我们还需要增加一个项目变量，因为不同的项目之间的差异也是需要考虑的。本例考虑每一个被试在每个条件下有12次观测。
#添加项目变量item for(i in 1:6){ item }oats
构建模型
建立混合线性模型：使用lme4包中的lmer函数
lmer(data = , formula = y ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor))
参数说明： data ：数据集 formula ：R中的公式 y :因变量 Fixed_Factor ：固定因子或自变量 Random_intercept ：随机截距。不同个体因变量分布不同。例如：有的被试反应程度比较快，有的比较慢。模型中一般是使用1表示，如果说不需药随机截距的话将1改为-1即可。 Random_Slope ：随机斜率。不同个体自变量与因变量关系不同。例如：有的被试对噪声有敏感，同时有的会不敏感。 Random_Factor ：随机因子。例如：两种实验中的被试、项目等等。
这里出现报错了，也就是随机效应的个数(144)大于样本数(72)。所以在建立模型过程，将过多的随机效应增加到模型中，当随机效应的个数大于样本数时会报错，这时需适当的调整、取消一些随机效应或者增加样本数。
题外话
表示即考虑固定因子的效应，也考虑其交互效应
+ 表示只考虑固定因子的效应，
：表示只考虑其交互效应
如果说对formula这个参数的设置不清楚，可以去看线性回归部分。点击转换
为了能继续下去，就需要对随机效应进行虚拟变量编码，有选择性的对随机效应选取。注：模型会随着加入的随机效应的不同而不同。这里仅为了展示所以选择随机效应时比较随意，真实中需重复试验考察。
展示一下编码后的数据(由于数据太多了，仅展示数据结构)
由原本的4列变成了17列，增加13列。包括item(项目变量)以及从(intercept)开始后面的12项随机效应的。
构建虚拟变量具体代码如下：
mm head(mm)data View(data)
同时为了方便，就建立一个零模型(不考虑协方差的作用)，将被试和项目前面的竖线改为两条即可。以下为更改后的模型。
oatslmer = lmer(yields ~ varietiesnitrogen + (1+varieties1+varieties2+nitrogen1+nitrogen2+nitrogen3+varieties1:nitrogen1+varieties2:nitrogen1+varieties1:nitrogen2+varieties2:nitrogen2+varieties1:nitrogen3+varieties2:nitrogen3 || blocks)+ (1+varieties1+varieties2+nitrogen1+nitrogen2+nitrogen3+varieties1:nitrogen1+varieties2:nitrogen1+varieties1:nitrogen2+varieties2:nitrogen2+varieties1:nitrogen3+varieties2:nitrogen3 || item), data = data)
更改后的模型，尽管模型没有报错，但是还是出现警告了。这个警告的意思就是模型中出现了畸形协方差矩阵。一般情况下零是不会出现畸形协方差矩阵的，因为零模型不考虑协方差的作用。尽管有警告也没影响。
畸形协方差是什么意思呢？我列举一个简单的例子(如下所示的模型)：
该模型也出现了畸形协方差矩阵。通过去观察随机效应可以得到为什么会存在警告。
在summary详细信息里面的randomeffects(随机效应)，最右边的corr就表示协方差，可以看到有一些相关性是非常强的，有的还等于1了，这就是协方差矩阵了。
如何判断模 型是否存在畸形协方差？
通过函数isSingular()检查，若返回为TRUE，则说明模型存在畸形协方差。
isSingular(oatslmer)[1] TRUE
注：模型有时候可能还会存在一种警告：模型会出现无法收敛的情况，也是通过调整或者取消一部分的随机效应。
回归正题
那么回到我们的零模型当中， 当出 现 这些问题的时候我们 要 如何去解决或者优化模型呢？
可通过主成分分析，确认有效的随机效应的成分数量，也就是剔除哪些对模型没有影响的随机效应成分。
在item随机因子上的主成分效应，第二行表示解释一个随机截距和十一个随机斜率的效应的方差，从第4个主成分开始能解释的比例为0了，因此在item随机因子上有9个随机效应时多余的。虽然后面的4和5仍有贡献，但是很小就不考虑进来了，选择剔除了，但是具体哪些是多余的，目前还不清楚。同理在blocks随机因子上，第9个主成分开始能够解释的比例为0了，所以在blocks随机因子上有4个随机效应是多余的。在本次建模剔除item随机因子上的9个随机效应和blocks随机因子上的4个随机效应。
怎么确定需要剔除的随机效应呢？
通过观察随机效应的标准差有两种方式：
一：通过查看原始模型中的Random effects中的StdDev
二：通过函数VarCorr()
两者的结果是一样的。
通过前面的主成分分析在item随机因子上剔除9个随机效应，也就是剔除标准差比较小的随机效应。同理在bolcks随机因子上剔除4个随机效应。
重新建模以后发现仍存在畸形协方差，再次采用主成分优化，经过N次以后得到最终的模型即：仅考虑blocks随机因子下的随机截距，item随机因子下的品种2和施肥2的交互与品种2和施肥3的交互随机效应。具体模型如下：
oatslmer1 = lmer(yields ~ varietiesnitrogen + (1| blocks)+ (-1+varieties2:nitrogen2+varieties2:nitrogen3 || item), data = data)
尽管模型的没有报错和警告，但是模型的效果并不是很好吧！(可以尝试将其他随机效应放入模型中)
模型解读
固定效应是通过summary(model)查看，其有四个部分：
一：Scaled residuals：包含标准差的最小值、中位数、四分位数。
二：Random effects：随机效应部分
三：Fixed effects：固定效应部分
四：Correlation of Fixed Effects：固定效应的协方差。
固定效应的结果如上图所示，这里是将varieties=Goldenrain和nitrogen=00cwt作为参照。可以看出nitrogen1、nitrogen2、nitrogen3(02cwt/04cwt/06cwt)其在参数估计(Estimate)要显著高于00cwt)。在品种varieties上和两者的交互作用上却没有显著差异。其余部分可自行去跑一下，这里由于篇幅问题就不展示了。
主效应和交互作用是通过anova(model)查看：
> anova(oatslmer1)Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) varieties 15306 7653 2 51918 35552 003572 nitrogen 181460 60487 3 50194 280982 8313e-11 varieties:nitrogen 2775 462 6 38619 02148 096984 ---Signif codes: 0 ‘’ 0001 ‘’ 001 ‘’ 005 ‘’ 01 ‘ ’ 1
结果如上，可看出燕麦的品种和施肥处理上的主效应显著，但其交互作用都不显著。
最后我们通过方差分析判断两个模型是否存在差异
方差分析的P值等于08，也就是模型间不存在显著性差异。即两个模型是一样的。
事后检验
模型的主效应显著后需要进行事后比较
由于我们模型的交互作用并不显著仅是主效应显著，那么对主效应进行多重比较。首先对品种(varieties)做多重比较。结果如下：
> emmeans(oatslmer1,pairwise~varieties,adjust="none")$emmeans varieties emmean SE df lowerCL upperCL Goldenrain 1044 703 653 875 121 Marvellous 1094 700 647 926 126 Victory 974 722 693 803 115Results are averaged over the levels of: nitrogen Degrees-of-freedom method: kenward-roger Confidence level used: 095 $contrasts contrast estimate SE df tratio pvalue Goldenrain - Marvellous -498 450 496 -1106 02742 Goldenrain - Victory 698 477 549 1462 01495 Marvellous - Victory 1196 467 502 2562 00134 Results are averaged over the levels of: nitrogen Degrees-of-freedom method: kenward-roger
contrasts返回的就是多重比较的结果。可以看到只有品种(marvellous)与品种(victory)是存在显著差异的。也可以通过emmeans返回的结果中可以看到品种(marvellous)燕麦的产量均值是1094而品种(victory)的产量均值为974。
我们再来看看施肥处理(nitrogen)的多重比较。结果如下：
> emmeans(oatslmer1,pairwise~nitrogen,adjust="none")$emmeans nitrogen emmean SE df lowerCL upperCL 00cwt 793 726 737 624 963 02cwt 991 723 731 821 1160 04cwt 1133 735 761 962 1304 06cwt 1233 745 781 1060 1405Results are averaged over the levels of: varieties Degrees-of-freedom method: kenward-roger Confidence level used: 095 $contrasts contrast estimate SE df tratio pvalue 00cwt - 02cwt -197 517 491 -3819 00004 00cwt - 04cwt -339 526 506 -6452 <0001> 00cwt - 06cwt -439 529 504 -8314 <0001> 02cwt - 04cwt -142 530 473 -2680 00101 02cwt - 06cwt -242 543 457 -4461 00001 04cwt - 06cwt -100 554 543 -1807 00763 Results are averaged over the levels of: varieties Degrees-of-freedom method: kenward-roger
同理contrasts返回的是施肥处理的多重比较结果。从上面可以看到除了04cwt与06cwt不存在显著差异外，其余的均存在显著性差异。
总结
一般来说做事后检验时，首先看交互作用，交互作用显著的话，就需要进行简单分析，简单分析中主效应显著了，再做多重比较。若交互作用不显著，那就直接对主效应做多重比较。也就是说，当你交互作用显著时，那么单纯做主效应分析的意义不大了，总之，交互作用的显著性优先考虑。(虽然该模型的交互作用不显著，但是还是走一遍流程吧！)
例子：
> #简单效应分析:比较不同品种上施肥处理的主效应> joint_tests(oatslmer1,by="varieties")varieties = Goldenrain: model term df1 df2 Fratio pvalue nitrogen 3 4031 9431 00001 varieties = Marvellous: model term df1 df2 Fratio pvalue nitrogen 3 4220 7063 00006 varieties = Victory: model term df1 df2 Fratio pvalue nitrogen 3 1494 8103 00019
结果：varieties = Goldenrain时，施肥处理的主效应是显著的；同时其他两种的品种在施肥处理上的主效应是显著的。
那么我的主效应显著以后就需要做多重比较了。
> emmeans(oatslmer1,pairwise~nitrogen|varieties,adjust="none")$emmeansvarieties = Goldenrain: nitrogen emmean SE df lowerCL upperCL 00cwt 802 880 149 615 990 02cwt 982 876 148 795 1168 04cwt 1131 927 166 935 1327 06cwt 1262 919 164 1067 1456varieties = Marvellous: nitrogen emmean SE df lowerCL upperCL 00cwt 864 879 150 677 1051 02cwt 1086 876 148 899 1273 04cwt 1171 902 161 980 1362 06cwt 1255 916 164 1061 1449varieties = Victory: nitrogen emmean SE df lowerCL upperCL 00cwt 714 912 166 521 907 02cwt 905 926 170 710 1100 04cwt 1097 1003 140 881 1312 06cwt 1182 1055 130 954 1410Degrees-of-freedom method: kenward-roger Confidence level used: 095 $contrastsvarieties = Goldenrain: contrast estimate SE df tratio pvalue 00cwt - 02cwt -1794 857 403 -2093 00427 00cwt - 04cwt -3288 921 439 -3571 00009 00cwt - 06cwt -4596 918 434 -5007 <0001> 02cwt - 04cwt -1494 918 440 -1628 01107 02cwt - 06cwt -2801 917 424 -3055 00039 04cwt - 06cwt -1307 982 457 -1331 01898 varieties = Marvellous: contrast estimate SE df tratio pvalue 00cwt - 02cwt -2217 865 422 -2563 00140 00cwt - 04cwt -3069 897 468 -3420 00013 00cwt - 06cwt -3911 914 448 -4279 00001 02cwt - 04cwt -852 898 462 -0949 03476 02cwt - 06cwt -1694 909 458 -1864 00687 04cwt - 06cwt -843 955 511 -0882 03817 varieties = Victory: contrast estimate SE df tratio pvalue 00cwt - 02cwt -1910 949 545 -2012 00492 00cwt - 04cwt -3824 1035 164 -3693 00019 00cwt - 06cwt -4677 1074 149 -4354 00006 02cwt - 04cwt -1915 1075 122 -1780 00998 02cwt - 06cwt -2767 1098 131 -2521 00255 04cwt - 06cwt -852 1188 127 -0717 04861 Degrees-of-freedom method: kenward-roger
式子：pairwise~nitrogen|varieties，其中：竖线前是预测变量，竖线后调节变量，如果实在不理解的话，你就理解是一个分式，竖线前是分子，竖线后是分母。当你要比较分子间的差异时，是不是要保证分母是相同的，也就是说，我要比较不同的施肥处理在燕麦产量上的差异时，我要保证的是施肥处理应作用在同一个品种上才可以进行比较？(在相同的品种上，不同施肥处理的方式间的差异是否显著)
结果还是看contrasts部分。例如varieties(品种) = Goldenrain上，施肥处理为02cwt - 04cwt和04cwt - 06cwt之间是不存在显著性(p值小于005)差异的，其他同理。
反之，如果你的简单主效应分析的是：不同施肥处理上品种的主效应，若显著后，需要做多重比较(在相同的施肥处理方式上，不同品种间的差异是否显著)
注：多重比较有很多比较的方法如下：
> padjustmethods[1] "holm" "hochberg" "hommel" "bonferroni" "BH" "BY" "fdr" [8] "none"
拓展
一：混合线性模型不仅可以使用lme4包中lmer()函数，也可以使用nlme包中lme()函数，这两个包语法都差不多，相对来说lme4这个包的运算速度会快一点。还有一个ASReml-R包中asreml()函数。
二：如果因变量是分类变量，则需要用广义线性模型。
glmer(data = , formula = , family = ,)其中，data就是数据集，formula和上面是一样的，family和逻辑回归中是一样的。
代码
#释放内存rm(list=ls())gc()#加载数据：采用R的MASS包中的oats数据集library(MASS)data(oats)names(oats) = c('blocks', 'varieties', 'nitrogen', 'yields')str(oats)#虚拟变量编码方式library(lme4)contrasts(oats$varieties)contrtreatment(3)contrsum(3)contrpoly(3)contrasts(oats$nitrogen)table(oats$varieties)table(oats$nitrogen)DV function(n){ m 1) for(i in 1:n){ for(j in 1:n-1){ if(i-1 == j) m[i,j] 1-( else m[i,j] 1/n) } } return(m)} DV(2)#重置品种和施肥处理的虚拟变量编码contrasts(oats$varieties) 3)contrasts(oats$nitrogen) 4)#添加项目变量item for(i in 1:6){ item 1:}oats #创建虚拟变量mm View(mm)data View(data)str(oats)str(data)#建模oatslmer_full = lmer(yields ~ varietiesnitrogen + (1+varietiesnitrogen | blocks)+ (1+varietiesnitrogen | item), data = oats)#模型优化oatslmer = lmer(yields ~ varietiesnitrogen + (1+varieties1+varieties2+nitrogen1+nitrogen2+nitrogen3+varieties1:nitrogen1+varieties2:nitrogen1+varieties1:nitrogen2+varieties2:nitrogen2+varieties1:nitrogen3+varieties2:nitrogen3 || blocks)+ (1+varieties1+varieties2+nitrogen1+nitrogen2+nitrogen3+varieties1:nitrogen1+varieties2:nitrogen1+varieties1:nitrogen2+varieties2:nitrogen2+varieties1:nitrogen3+varieties2:nitrogen3 || item), data = data)summary(oatslmer)isSingular(oatslmer)summary(rePCA(oatslmer))VarCorr(oatslmer)#最终模型oatslmer1 = lmer(yields ~ varietiesnitrogen + (1| blocks)+ (-1+varieties2:nitrogen2+varieties2:nitrogen3 || item), data = data)summary(oatslmer1)anova(oatslmer1)anova(oatslmer,oatslmer1)#主效应事后比较library(emmeans)padjustmethodsemmeans(oatslmer1,pairwise~varieties,adjust="none")emmeans(oatslmer1,pairwise~nitrogen,adjust="none")#简单效应分析:比较不同品种上施肥处理的主效应joint_tests(oatslmer1,by="varieties")emmeans(oatslmer1,pairwise~nitrogen|varieties,adjust="none")
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）