这篇文章发表于 1772 天前，可能其部分内容已经发生变化，如有疑问可询问作者。

因为之前看到不少资料，但好多文章缺乏整体的框架意识，导致我老是会有错乱的感觉，希望这份东西至少能够满足我自己的需求——没有任何多余的步骤、过分的解释和数学证明，单纯介绍一部分基本概念在R中的利用方式并努力让读者建立起一定的框架。文中仅涉及部分核心函数的介绍，至于函数功能的实现，文中大多仅展示部分，当然可以有其他办法。

整篇文章就是由各种资料拼接而成，感谢这些师傅和所有为之付出过心血的人。

本文的最后将列出一部分个人觉得非常好的文章。

估计

这部分主要是理论基础，没啥好记的。不过bootstrap比较好用。

1
2
#回归参数的区间估计
confint(object, parm, level = 0.95, ...)

假设检验

参数检验与非参数检验区别：参数检验假定数据服从某分布（一般为正态分布），通过样本参数的估计量对总体参数进行检验；而非参数检验不需要假定总体分布形式，不涉及总体分布的参数，直接对数据的分布进行检验。

在检验之前，预先对分布检验，也就是正态性的检验非常重要。另外，对独立性的检验也是很重要的，这可能需要秩相关分析（见第3部分）来确定。使用独立样本非参数检验还是相关样本非参数检验需要通过自己的专业判断。

一些较常用的假设检验摘录于此，还有一些检验由于排版需要，放在了后面部分。

a. 参数检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#Bartlette检验是方差齐性检验，对于正态分布的样本，Bartlette检验极其灵敏，但是对于非正态分布的样本，检验非常不准确
bartlett.test(x, g, ...)

#levene检验同样是方差齐性检验，需要导入car包，levene检验既可用于正态分布的样本，也可用于非正态分布的样本，同时对比较的各组样本量可以相等或不等
leveneTest(y, group, center=median, ...)
#levene.test函数已被弃

#Z检验，一般用于大样本(即样本容量大于30)，平均值差异性检验BSDA包提供了函数z.test()
z.test(x, y = NULL, 
       alternative = "two.sided", 
       mu = 0, sigma.x = NULL, sigma.y = NULL, 
       conf.level = 0.95)

#t检验主要用于样本含量较小（例如n < 30），总体标准差σ未知的正态分布（所以建议先检验分布类型）。t检验是用t分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著
t.test(x, y = NULL,
       alternative = c("two.sided", "less", "greater"),
       mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95, ...)
#alternative即H1，greater对应区域1、2顺序 ，反过来要把alternative改为less；另外var.equal在方差不等的时候设为FALSE
#我以前的文章，可参看这里

#两总体方差比较：F检验（应用条件：所有分组数据正态分布，并且方差要齐）
var.test(x, y, ratio = 1,
         alternative = c("two.sided", "less", "greater"),
         conf.level = 0.95, ...)
#F检验常被用于检验线性约束条件，另外还有（计量经济学三大）检验同样可以实现——似然比检验（LR）、沃尔德检验（Wald）和拉格朗日乘子检验（LM）。三者的操作细节参见底部的参考文献

b. 独立样本非参数检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#二项分布检验，检验二项分布概率，需要导入stats包
#原假设：样本来自概率为指定数值的二项分布总体
binom.test(x, n, p = 0.5,
           alternative = c("two.sided", "less", "greater"),
           conf.level = 0.95)

#游程检验，检验样本的随机性，也可检验两个总体的分布是否相同，游程总数不宜太多太少，需要导入tseries包
#原假设：样本具有随机性
runs.test(x, alternative = c("two.sided", "less", "greater"))

#卡方检验（皮尔森拟合优度卡方检验，可作独立性检验），用于检验数据的一致性，需要导入stats包
#原假设：某样本总体分布与另一分布无显著差异
chisq.test(x, y = NULL, correct = TRUE,
           p = rep(1/length(x), length(x)), rescale.p = FALSE,
           simulate.p.value = FALSE, B = 2000)
#correct用于表示是否连续修正，simulate.p.value以及B与仿真有关；另外，该函数支持列联表检验

#Two-proportions test用来检验两组或多组二项分布的成功概率（比例）是否相等，或等于给定的值，需要导入stats包
#该函数利用了小样本单尾卡方检验
#原假设：样本间无明显差异
prop.test(x, n, p = NULL,
          alternative = c("two.sided", "less", "greater"),
          conf.level = 0.95, correct = TRUE)

#S-W检验，即Shapiro-Wilk检验（常作为正态性检验）,适用于小样本场合，需要导入stats包
#原假设：样本满足正态分布
shapiro.test(x) 
#x a numeric vector of data values. x可以是数值型向量，允许存在NA，但是非丢失数据需要在3-5000内。

#K-S检验，即Kolmogorov-Smirnov检验（常作为正态性检验）,适用于大样本场合，尤其是连续分布。当分布参数未知时, 结果不可靠, 使用 Monte-Carlo方法，需要导入stats包
#原假设：某样本总体分布与另一分布无显著差异
ks.test(x, y, ..., 
        alternative = c("two.sided", "less", "greater"), 
        exact = NULL)

#Wilcoxon秩和检验（也称Mann-Whitney U检验），功能上类似于均值比较检验，常在服从正态分布的假设并不确定时使用。注意，样本中存在大量“结”（相同的数值）时，此法会出问题。需要导入stats包
#原假设：两样本来自总体的分布无显著差异
wilcox.test(x, y = NULL,
            alternative = c("two.sided", "less", "greater"),
            mu = 0, paired = FALSE, exact = NULL, correct = TRUE,
            conf.int = FALSE, conf.level = 0.95, ...)
#paired逻辑变量说明是否成对样本，exact是否精确计算p值，在x,y均存在并且paired=FALSE时，该函数自动选择秩和检验，反之为符号秩检验
#Density, distribution function, quantile function and random generation for the distribution of the Wilcoxon rank sum statistic obtained from samples with size m and n, respectively.
#dwilcox(x, m, n, log = FALSE)
#pwilcox(q, m, n, lower.tail = TRUE, log.p = FALSE)
#qwilcox(p, m, n, lower.tail = TRUE, log.p = FALSE)
#rwilcox(nn, m, n)

#Moses极端反应检验，一般用于检验实验组与对照组之间的差异是否显著，需要导入DescTools包
#原假设：某样本总体分布与另一分布无显著差异
MosesTest(x, y, extreme = NULL, …)

#Wald-Wolfowitz游程检验，考察样本的随机性，需要导入randtests包
#原假设：样本具有随机性。注意，在SPSS中，该函数貌似经过处理，作为检验两样本是否具有相同分布，其备选假设是两样本在某方面存在差别
runs.test(x, alternative, threshold, pvalue, plot)
#plot=TRUE|FALSE

#中位数检验（Mood's Median Test），检验多个总体的中位数是否存在显著差异，需要导入onewaytests包
#原假设：k个总体的中位数相同
mood.test(y, group, na.rm = TRUE)

#Jonckheere-Terpstra检验，一般用在k个总体已经自然先验排序后，需要导入clinfun包
#原假设：两样本来自同一总体
jonckheere.test(x, g, alternative = c("two.sided", "increasing",
                 "decreasing"), nperm=NULL)

#另外，Kruskal-Wallis H检验，相当于Mann-Whitney U检验的扩展
#原假设：两样本来自同一总体
#具体方法在ANOVA非参部分提及

c. 相关样本非参数检验

这里主要是两样本的方法，多样本的方法写在了非参数ANOVA部分。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#McNemar变化显著性检验，针对名义级变量（nominal），两组样本观察数目相同并且观察值顺序不能随便改，注意样本是二分类数据，需要导入stats包
#原假设：两组样本分布没有显著差异
mcnemar.test(x, y = NULL, correct = TRUE)

#边际同质性（Marginal Homogeneity）检验，相当于分类数据的McNemar检验，针对顺序级变量（ordinal），需要导入DescTools包
#原假设：两组样本分布没有显著差异
StuartMaxwellTest(x, y = NULL)

#符号（Sign）检验，检验两组相关数值型数据的总体分布是否相同，针对数值型变量（scale），需要导入BSDA包
#原假设：两组样本分布没有显著差异
SIGN.test(x, y = NULL, md = 0, alternative = "two.sided",
  conf.level = 0.95, ...)

#Wilcoxon符号秩检验（Wilcoxon signed rank test），检验两组相关数值型数据的总体分布是否相同，针对数值型变量（scale），注意样本须独立，而且得是连续型变量，需要导入stats包
#原假设：两组样本分布没有显著差异
wilcox.test(x, y = NULL,
            alternative = c("two.sided", "less", "greater"),
            mu = 0, paired = FALSE, exact = NULL, correct = TRUE,
            conf.int = FALSE, conf.level = 0.95, ...)
#该函数在x,y均存在并且paired=FALSE时，该函数自动选择秩和检验，反之为符号秩检验

#Cochran's Q检验，此检验相当于二元响应的Friedman检验，是McNemar检验对k样本情况的扩展,适用于名义级变量的检验，需要导入RVAideMemoire包
#原假设：多个相关二分变量具有相同的平均值
cochran.qtest(formula, data, alpha = 0.05, p.method = "fdr")

#Kendall W协同系数检验，检验多个相关总体变化趋势是否一致，适用于顺序型或数值型变量，需要导入irr包
#参数W（范围0~1）接近1说明协同性好
kendall(ratings, correct = FALSE)
#correct意为是否修正

#另外，Friedman卡方检验适用于顺序型变量或是（可转换为定距的）数值型变量，操作在后面ANOVA非参分析中提及
#原假设：所有相关样本总体分布相同|不同机构的评价标准是不相关的

相关性分析

a. 简单线性相关与秩相关

1
2
3
4
5
6
7
8
#对一般情况默认数据服从正态分布的，用pearson相关系数（连续性变量才可采用）；当资料不服从双变量正态分布或总体分布未知，或原始数据用等级表示时，宜用spearman或kendall
cor.test(x, y,
         alternative = c("two.sided", "less", "greater"),
         method = c("pearson", "kendall", "spearman"),
         exact = NULL, conf.level = 0.95, continuity = FALSE, ...)
#pearson's r（一种参数方法）对异常值敏感；
#spearman's rho适用于连续与离散型变量，对异常值不敏感，另外样本量过小时，其产生的p值非常不准确；
#kendall's tau定性判断，对异常值不敏感

b. 偏相关

存在与两个变量都有相关性的混淆变量时，相关系数往往结果具有误导性，于是有了偏相关。

偏相关与半偏相关区别：

The semipartial (or part) correlation statistic is similar to the partial correlation statistic. Both compare variations of two variables after certain factors are controlled for, but to calculate the semipartial correlation one holds the third variable constant for either X or Y but not both, whereas for the partial correlation one holds the third variable constant for both. The semipartial correlation compares the unique variation of one variable (having removed variation associated with the Z variable(s)), with the unfiltered variation of the other, while the partial correlation compares the unique variation of one variable to the unique variation of the other. 

1
2
3
4
5
6
7
#partial correlations 偏相关，需要导入ppcor包
pcor.test(x, y, z, method = c("pearson", "kendall", "spearman"))
#在给定z条件下，研究x~y的偏相关性

#Semi-partial correlation 半偏相关，需要导入ppcor包
spcor.test(x, y, z, method = c("pearson", "kendall", "spearman"))
#在某x,y变量（不妨为y）移除z条件后，研究x~y的偏相关性

c. 典型相关

简单来讲，就是将一堆要素综合之后化为两个矩阵，再求两个矩阵之间的关系。

后面的d和e部分其实是在这个基础上的延伸应用。其中的mantel's test在d部分中会继续探讨。

1
2
3
4
5
6
7
8
9
10
11
12
#Compute the canonical correlations between two data matrices
cancor(x, y, xcenter = TRUE, ycenter = TRUE)

#Mantel test是对两个矩阵相关关系的检验，其原假设是两个矩阵间没有相关关系；mantel.partial类似于偏相关，都需要导入vegan包
mantel(xdis, ydis, method="pearson", permutations=999, strata = NULL,
    na.rm = FALSE, parallel = getOption("mc.cores"))
mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 999, 
    strata = NULL, na.rm = FALSE, parallel = getOption("mc.cores"))

#The Knox test，R源代码参看https://github.com/jimhester/surveillance/blob/master/R/knox.R

#Jacquez’s test，R源代码参看https://www.niph.go.jp/soshiki/gijutsu/download/Rfunctions/func/Jacquez.test

另外，典型相关分析简称CCA（canonical correlation analysis），它的简称和后面会提到的典型对应分析（canonical correspondence analysis）相同。但是后者的目的是识别并量化两个数据集之间的相关关系，关注的是两数据集各自的线性组合的相关关系；典型相关分析测定的是两个数据集（多个列向量）之间的相关程度，整个过程中试图把两高维数据的相关关系压缩到少数几个典型变量对之间。

d. 时空自相关

这主要用两个矩阵，一个是时间的，另一个是空间的，再用mantel's test来检验这两个矩阵的关系。以下是mantel’s test的大致原理：

These are the two matrices which the function will be testing for a correlation. The test consists of calculating the correlation of the entries in the matrices, then permuting the matrices and calculating the same test statistic under each permutation and comparing the original test statistic to the distribution of test statistics from the permutations to generate a p-value. The number of permutations defines the precision with which the p-value can be calculated.

1
2
3
#Mantel’s test 作为空间自相关的检验也可（但在这里主要是考虑空间距离矩阵G和时间距离矩阵），需要导入ade4包
mantel.rtest(m1, m2, nrepet = 99, ...)
#其中m1，m2是矩阵（比如m1是时间距离矩阵，m2是空间距离矩阵，如此构造dist(cbind(LON,LAT))即可），而nrepet是置换次数，决定了精度

e. 空间自相关

自己补了一篇Moran's I操作方面的简易指北文章，在这里。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#Moran's I（全局） 以及 Local Moran’s I - Moran’s Ii（局部），需要导入spdep包
moran.test(x, listw, randomisation=TRUE, zero.policy=NULL,
 alternative="greater", rank = FALSE, na.action=na.fail, spChk=NULL,
 adjust.n=TRUE, drop.EI2=FALSE)
#这里的空间权重矩阵listw可以用nb2listw由nb对象（neighbours list object）转为listw对象（listw是在空间回归命令中可以直接使用的一种空间矩阵形式）

#制作Moran's I的散点图
moran.plot(x, listw, zero.policy=NULL, spChk=NULL, labels=NULL,
 xlab=NULL, ylab=NULL, quiet=NULL, ...)

#Getis-Ord General G 方法，需要导入spdep包
globalG.test(x, listw, zero.policy=NULL, alternative="greater",
 spChk=NULL, adjust.n=TRUE, B1correct=TRUE, adjust.x=TRUE, Arc_all_x=FALSE)

#join-count方法，二值或是类型条件下使用，需要导入spatial包
joincount.multi(fx, listw, zero.policy = FALSE,
 spChk = NULL, adjust.n=TRUE)
#也有函数joincount.test，但推荐multi

方差分析ANOVA

ANOVA表格就不放了。[doge] 这部分本文涉及不多，建议多看看链接。

原假设：模型间不存在显著差异。

a. 普通ANOVA

前提条件：各因素水平下的数据服从正态分布并且满足方差齐性。这些检验方式详见文章第2部分——假设检验。

1
2
3
4
#方差分析是从均值比较检验（单因素，双水平）扩展而来，有单因素方差分析、多因素方差分析（交互作用不明显）、多因素方差分析（含交互作用）、协方差分析，需要导入stats包
aov(formula, data = NULL, projections = FALSE, qr = TRUE,
    contrasts = NULL, ...)
#具体操作可参见http://biotrainee.com/jmzeng/markdown/ANOVA.html

b. 非参数ANOVA

在ANOVA所需的前提条件不满足时酌情使用替代方案，因为可能结果可能出现一定的问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
#Kruskal-Wallis检验，非参数方法，检验多个独立总体平均秩是否存在差异，可作为单因素方差分析的替代，需要导入stats包
#统计量W（范围0~1）接近1，就说明一致程度较高
#原假设：不同变量造成的效果无差异
kruskal.test(x, g, …)

#Friedman卡方检验，非参数方法，可作为双因素方差分析 （无交互作用）的替代，需要导入stats包
#原假设：所有相关样本总体分布相同
friedman.test(y, groups, blocks, …)
#It’s recommended when the normality assumptions of the one-way repeated measures ANOVA test is not met or when the dependent variable is measured on an ordinal scale.

#Scheirer-Ray-Hare检验，可作为多因素方差分析 （含交互作用）的替代
scheirerRayHare(formula = NULL,data = NULL,y = NULL,x1 = NULL,x2 = NULL,
                tie.correct = TRUE,ss = TRUE,verbose = TRUE)

回归分析

回归分析目的是根据已知的自变量来估计或者预测因变量的情况 。说实话以前觉得这里的数学推导并不难，但是前提假设上有各种限制与坑点，这部分并不简单（要是以后有空有心情就写一篇完整的文章阐述一下），所以这里我要把一些东西写得稍微详细一点。

由于限制过多，可先尝试回归后再回过头检验，但按逻辑，还是将这些限制与相关检验写在pre部分。

pre0. 由假设产生的限制

在线性回归模型中，残差（即随机误差）应满足白噪声假设（White Noise Condition）：1、残差（无自相关性）独立同分布；2、残差和自变量不相关；3、残差均值为0，方差为常数。

ps：理论上，白噪声假设不要求随机变量服从正态分布，而可以是任意分布。但基于中心极限定理，假设残差服从正态分布是一个合理的近似。

由此，需要对假设条件进行检验，一般分为这些：

  • 自相关性（也称序列相关性）：若自相关性较高，在其他条件好的情况下，参数估计依旧无效——因为方差存在偏误且均值非无偏估计。检验方法有Breusch－Pagan test、White test，处理方法主要使用广义差分法。
  • 异方差性：即残差方差不为常数。适合回归的应不存在异方差性。其不影响最小二乘估计值的无偏性，但影响有效性和显著性。检验方法有White test，处理方法有加权最小二乘法（weighted least squares, WLS）。
  • 多重共线性：自变量间存在较大的相关性，会导致t检验不可靠，很可能得出非常不合理的结果。主要判别依据有VIF、容忍度、CI(Conditional Index)等。可以通过剔除不必要自变量、主成分回归、岭回归以及改变自变量或者函数形式这些方法尽可能消除其影响。
  • 内生性：源于OLS（ordinary least squares 普通最小二乘法）其中的解释变量与随机扰动间存在联系。在没有内生性问题的情况下，OLS更为合适。检验方法有Hausman test。后面对这一部分不予讨论。

对一些违背假设的模型依旧可以通过上述针对性方法处理后来回归分析。

pre1. 自变量优选

前进法、后退法与逐步回归法的操作方法可参见这里以及这里。

pre2. 对利用回归前提的检验

检验自相关性：

1
2
3
#Durbin-Watson test，函数名可被缩写成dwt，需要导入car包
durbinWatsonTest(model, max.lag=1, ...)
#D.W统计量在2左右说明残差是服从正态分布的，若偏离2太远，靠近0或4，那么就存在自相关性。其p值计算时原假设是不存在自相关。

检验异方差性：

1
2
3
4
5
6
7
8
#Breusch－Pagan test，需要导入lmtest包
bptest(formula, varformula = NULL, studentize = TRUE, data = list())

#white test其实可看作是广义bptest的一种特殊形式，故能够通过修改参数后进行bptest来实现，其原假设是不存在异方差性
#比如回归为fm <- lm(y ~ x + z, data = foo)，那么则应写成bptest(fm, ~ x * z + I(x^2) + I(z^2), data = foo) 
#也可以通过导入tseries包实现
white.test(x, y, qstar = 2, q = 10, range = 4,
            type = c("Chisq","F"), scale = TRUE, ...)

检验多重共线性：

1
2
3
#方差膨胀因子VIF，是容忍度（tolerance）的倒数，需要导入car包
vif(mod, ...)
#如果全部小于10（严格是5），则说明模型没有多重共线性问题，模型构建良好；反之若VIF大于10说明模型构建较差。

a. 简单线性回归/多元线性回归

1
2
3
4
5
6
#线性回归方程的原假设：因变量与自变量之间无显著关系
#线性回归系数的原假设：beta[i]=0，也就是该自变量与因变量无显著关系
lm(formula, data, subset, weights, na.action,
   method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
   singular.ok = TRUE, contrasts = NULL, offset, ...)
#简单线性回归y~x，多元线性回归y~x1+x2+...，该函数一般结合函数abline与summary

b. 广义线性回归GLM

这里包含了各种含定性（以及定序）变量的回归。其估计采用最大似然法 ，参数显著性检验采用Wald test (Z 统计)，模型显著性使用似然比检验 。

1
2
3
4
5
6
7
8
#广义线性回归，需要导入stats包
glm(formula, family = gaussian, data, weights, subset,
    na.action, start = NULL, etastart, mustart, offset,
    control = list(...), model = TRUE, method = "glm.fit",
    x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)
#参数family可选[quasi]binomial(link="logit")，[quasi]poisson(link = "log"),gaussian(link = "identity"),Gamma(link = "inverse"),inverse.gaussian(link = "1/mu^2"),quasi(link = "identity", variance = "constant")
#这里使用binomial实现logistic回归；使用poisson实现Poisson回归
#因变量是二值型变量(Binary), 类别/层次型变量(Categorical)和次序变量(Ordinal)时，宜采用Logistic回归；因变量是计数变量(count)时, 宜使用Poisson回归

c. 非线性回归

非线性回归是回归函数关于未知回归系数具有非线性结构的回归，可分为可线性化问题和不可线性化问题。常用的处理方法有回归函数的线性迭代法、分段回归法、非线性最小二乘法（NLS）等。

1
2
3
4
#非线性最小二乘法使用线性函数来逼近非线性函数，并且通过不断迭代这个过程来得到参数的最优解，需要导入stats包
nls(formula, data, start, control, algorithm,
    trace, subset, weights, na.action, model,
    lower, upper, ...)

d. 稳健回归

稳健回归(robust regression)是将稳健估计方法用于回归模型，以拟合大部分数据存在的结构，同时可识别出潜在可能的离群点、强影响点或与模型假设相偏离的结构。有这样的一些方法（部分是非参数方法）：Brown-Mood估计、Theil-Sen回归（或Theil估计） 、Siegel估计 、S估计、小中位数二乘（LMS）、小截尾二乘（Least Trimmed Squares, LTS）、线性分位数回归 、Huber回归 、随机采样一致性算法（RANdomSAmpleConsensus, RANSAC）。

一般的稳健性排序： OLS< Theil < Siegel < (S, LTS, LMS)

e. 分位数回归

当不满足随机误差的“正态独立同分布”，或数据波动较大时，可降格采用分位数回归（Quantile regression），属于非参数统计。

1
2
3
4
#分位数回归，需要安装quantreg包
rq(formula, tau=.5, data, subset, weights, na.action, method="br", 
   model = TRUE, contrasts, ...)
#使用方式可参见https://www.bbsmax.com/A/MAzAoxy59p/

f. 线性混合效应模型

混合效应分为固定(fixed)效应和随机(random)效应两部分。固定效应是指经典线性回归的部分, 随机效用是指分组造成的差异部分（随机这个名称与数学随机性没有太大关系） 。在许多情况下，同一变量可能被视为随机效应或固定效应（有时甚至是同时发生！），因此请始终参考您的问题和假设以相应地构建模型。

对于嵌套随机效应，该因子仅出现在另一个因子的特定级别内； 对于交叉效应，给定因子出现在另一个因子的一个以上层次中。 或者，您还可以记住，如果您的随机效果没有嵌套，那么它们就会交叉！

混合模型使我们可以节省自由度。

1
2
3
4
5
6
7
8
9
10
11
#线性混合效应模型，需要安装lme4包
lmer(formula, data = NULL, REML = TRUE, control = lmerControl(),
     start = NULL, verbose = 0L, subset, weights, na.action,
     offset, contrasts = NULL, devFunOnly = FALSE, ...)
#强烈推荐https://ourcodingclub.github.io/tutorials/mixed-models/
#使用以上网站中的数据，另注：
#dragons <- within(dragons.data, sample <- factor( mountain : site))
#mixed.lmer <- lmer(score ~ length + (1 | mountain) + (1 | sample), data = dragons)
#与
#fit3 <- lmer(score ~ 1 + length + (1 | mountain / site), data=dragons.data) ##常数项取0或-1时表示零截距回归，不写或取1表示此处不强制使截距为零
#所得结果相同

g. 地理加权回归

1
2
3
4
5
6
7
8
9
10
11
#地理加权回归GWR，需要导入spgwr包
#可以选出带宽看看
gwr.sel(formula, data=list(), coords, adapt=FALSE, gweight=gwr.Gauss,
 method = "cv", verbose = TRUE, longlat=NULL, RMSE=FALSE, weights,
 tol=.Machine$double.eps^0.25, show.error.messages = FALSE)
#这里的coords一般是cbind(X,Y)
#回归
gwr(formula, data=list(), coords, bandwidth, gweight=gwr.Gauss, 
	adapt=NULL, hatmatrix = FALSE, fit.points, longlat=NULL, 
	se.fit=FALSE, weights, cl=NULL, predictions = FALSE, 
        fittedGWRobject = NULL, se.fit.CCT = TRUE)

h. 空间自回归

空间自回归模型(Spatial autoregressive model)是指因变量线性依赖于自身邻居的作用、 其它解释变量的共同作用和随机成分。可结合使用AIC来选取较好的模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#以下操作需要导入spatialreg包
#lag回归模型，y=ρWy+Xβ+ε
lagsarlm(formula, data = list(), listw, na.action, Durbin, type,
 method="eigen", quiet=NULL, zero.policy=NULL, interval=NULL,
 tol.solve=.Machine$double.eps, trs=NULL, control=list())
#error回归模型，y=Xβ+u,u=λWu+ε
errorsarlm(formula, data=list(), listw, na.action, weights=NULL,
 Durbin, etype, method="eigen", quiet=NULL, zero.policy=NULL,
 interval = NULL, tol.solve=.Machine$double.eps, trs=NULL,  control=list())
#sac回归模型，y=ρW1y+Xβ+u,u=λW2u+ε
sacsarlm(formula, data = list(), listw, listw2 = NULL, na.action, Durbin, type,
 method="eigen", quiet=NULL, zero.policy=NULL, tol.solve=.Machine$double.eps,
 llprof=NULL, interval1=NULL, interval2=NULL, trs1=NULL, trs2=NULL,
 control = list())

i. Spatial Panel Data Models

是我飘了。这部分太难了。不知有无大哥愿意贡献些？

降维处理

a. 聚类分析

主要的层次聚类方法有以下几种——基于连接的层次聚类：单连接聚合聚类（单联动）、完全连接聚合聚类（全联动）；平均聚合聚类：使用算术平均的非权重成对组法（UPGMA）、使用质心的非权重成对组法（UPGMC）、使用算术平均的权重成对组法（WPGMA）、使用质心的权重成对组法（WPGMC）；Ward最小方差聚类。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#hclust聚类分析，需要导入stats包
hclust(d, method = "complete", members = NULL)
#其中method可选"ward.D", "ward.D2", "single" (= 单联动), "complete" (= 全联动), "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC)
#通常d是距离矩阵，可由dist函数或是vegan包中的vegist()转换数据而来
#作图，这里的函数算是被这个包劫持了？
plot(x, labels = NULL, hang = 0.1, check = TRUE,
     axes = TRUE, frame.plot = FALSE, ann = TRUE,
     main = "Cluster Dendrogram",
     sub = NULL, xlab = NULL, ylab = "Height", ...)

#kmeans聚类分析，只能处理多凸边形，需要导入stats包
kmeans(x, centers, iter.max = 10, nstart = 1,
       algorithm = c("Hartigan-Wong", "Lloyd", "Forgy",
                     "MacQueen"), trace=FALSE)
#centers为中心个数，nstart为重复次数

#Density-Based Spatial Clustering of Applications with Noise (dbscan)，可对非凸多边形进行处理，需要导入
dbscan(data, eps, MinPts = 5, scale = FALSE, method = c("hybrid", "raw",
    "dist"), seeds = TRUE, showplot = FALSE, countmode = NULL)
#eps是半径，MinPts最小点数默认为5
#另外，eps可通过dbscan::kNNdistplot来寻找，需要导入dbscan包
kNNdistplot(x, k, all = FALSE, ...)
#其中k即为MinPoints

b. 主成分分析PCA/因子分析FA

主成分分析是因子分析的一种。首先应该注意到因子分析的应用条件：

  • 样本量不能太小，至少为变量数的5倍。
  • 各变量间应该具有相关性，如彼此独立，则无法提取公因子，最好偏相关性也比较高。
  • 因子分析中各公因子应该具有实际意义。

确认待分析的原始变量之间是否存在较强的相关关系：

1
2
3
4
5
6
7
8
9
#KMO检验，用于考察变量间的偏相关性，取值0~1之间，需要导入psych包
KMO(r)
#其中r为 A correlation matrix or a data matrix (correlations will be found)
#KMO统计量越接近1，变量间的偏相关性越强，因子分析效果越好。一般统计量在0.7以上为适应做因子分析。＜0.5则不适宜做因子分析

#Bartlett球性检验，可参见前面参数检验部分内容，也可参考此函数，需要导入psych包
cortest.bartlett(R, n = NULL,diag=TRUE)
#其中R为 A correlation matrix.
#其原假设是原始变量间不存在相关性

实现PCA的简单方法，至于FA，可参见文末链接。

1
2
3
4
5
6
7
8
9
10
11
12
13
#Cattell碎石检验绘制，以确定主成分的构成，根据Kaiser-Harris准则建议保留特征值大于1的主成分，需要导入psych包
fa.parallel(x,n.obs=NULL,fm="minres",fa="both",nfactors=1, 
            main="Parallel Analysis Scree Plots", 
            n.iter=20,error.bars=FALSE,se.bars=FALSE,SMC=FALSE,ylabel=NULL,
            show.legend=TRUE,sim=TRUE,quant=.95,cor="cor",
            use="pairwise",plot=TRUE,correct=.5)

#PCA分析，需要导入psych包
principal(r, nfactors = 1, residuals = FALSE,rotate="varimax",n.obs=NA, 
          covar=FALSE,scores=TRUE,missing=FALSE,impute="median",
          oblique.scores=TRUE,method="regression",use="pairwise",
          cor="cor",correct=.5,weight=NULL,...)
#"none", "varimax", "quartimax", "promax", "oblimin", "simplimax", and "cluster" are possible rotations/transformations of the solution.

c. 对应分析CA

CA (Correspondence Analysis) 是R型与Q型FA的结合，是一种描述性统计，对极端值敏感。主要有简单对应分析（两个分类变量）、多重对应分析（多于两个分类变量）、数值变量对应分析、均值对应分析。

它虽然可以揭示变量间的联系，但不能说明两个变量之间的联系是否显著，因而在做对应分析前，可以用其他方法检验变量的相关性。原始数据的无量纲化处理，意味着处理时应注意各变量量纲一致。

1
2
3
4
5
6
7
8
9
#CA的一种流程方法，需要MASS包
corresp(x, ...)
#该函数作出了相关（得分）矩阵

#利用该函数对得出的矩阵作图，需要stats包
biplot(x, y, var.axes = TRUE, col, cex = rep(par("cex"), 2),
       xlabs = NULL, ylabs = NULL, expand = 1,
       xlim  = NULL, ylim  = NULL, arrow.len = 0.1,
       main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ...)

CA有时会产生“弓形效应”，对排序的精度产生影响，可通过去趋势对应分析（Detrended Correspondence Analysis，DCA）来解决，参见这里。

另外，典范对应分析 CCA 基于对应分析（CA）发展而来的一种排序方法，将对应分析与多元回归分析相结合，又称多元直接梯度分析，参见这里。

extra

这里记录一些常用而不适合放在前面部分的知识。

a. 不均衡系数

1
2
3
4
5
6
7
8
#需要导入ineq包
Gini(x, corr = FALSE, na.rm = TRUE) #结合Lc函数可作出洛伦茨曲线
RS(x, na.rm = TRUE)
Atkinson(x, parameter = 0.5, na.rm = TRUE)
Theil(x, parameter = 0, na.rm = TRUE) #如果Theil中的参数为0，则计算Theil的熵度量，对于其他每个值，计算Theil的第二度量
Kolm(x, parameter = 1, na.rm = TRUE)
var.coeff(x, square = FALSE, na.rm = TRUE)
entropy(x, parameter = 0.5, na.rm = TRUE) #可使用此函数代替Theil

b. plot函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
plot(x, y = NULL, type = "p",  xlim = NULL, ylim = NULL,
     log = "", main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
     ann = par("ann"), axes = TRUE, frame.plot = axes,
     panel.first = NULL, panel.last = NULL, asp = NA,
     xgap.axis = NA, ygap.axis = NA,
     ...)
#pch（plotting characters）指定绘制点时使用的符号;
#cex指定符号的大小。cex是一个数值，表示绘图符号相对于默认大小的缩放倍数（默认大小为1）1.5表示放大为默认值的1.5倍；
#lty（line types）指定线条类型 (0=blank, 1=solid (default), 2=dashed, 3=dotted, 4=dotdash, 5=longdash, 6=twodash) ；
#lwd（line widths）指定线条宽度；
#col（colors for lines and points）指定颜色；
#xlab指定y轴标题；
#type："p" for points,"l" for lines,"b" for both,"c" for the lines part alone of "b","o" for both ‘overplotted’,"h" for ‘histogram’ like (or ‘high-density’) vertical lines,"s" for stair steps,"S" for other steps, see ‘Details’ below,"n" for no plotting；
#asp（aspect ratio y/x）注：在使用经纬图时推荐大于1；
#las（style of axis labels）指定坐标刻度值文字方向，las=0表示文字方向与坐标轴平行，1表示始终为水平方向，2表示与坐标轴垂直，3表示终为垂直方向；
#axes在FALSE情况下隐藏坐标轴，单个的话使用xaxt,yaxt。
#xlim限制区间，如果前者比后者大，那么就会反转坐标轴
#ann指定题目、坐标轴等是否显示

*相关文章资料*

大综合：

https://mgimond.github.io/Spatial/index.html

https://wangcc.me/LSHTMlearningnote/

https://bookdown.org/yufree/datadown/

https://bookdown.org/yufree/Notes/

https://bookdown.org/hezhijian/book/

https://rspatial.org/raster/index.html

https://bookdown.org/tpemartin/econometric_analysis/

https://bookdown.org/xiangyun/masr/ （好像被作者删了大部分）

http://www.tinysoft.com.cn/tsdn/helpdoc/display.tsl?id=12856

https://otexts.com/fppcn/

https://cran.r-project.org/doc/contrib/Xu-Statistics_and_R.pdf

https://www.ibm.com/support/knowledgecenter/zh/SSLVMB_sub/statistics_mainhelp_ddita/spss/base/core_container.html （SPSS手册）

R的基本技能：

https://openbiox.github.io/Cookbook-for-R-Chinese/

http://www.cookbook-r.com/

https://www.jianshu.com/p/baa6b9da31ba （有关R的plot函数）

https://ggplot2.tidyverse.org/reference/index.html （有关R的ggplot函数）

http://blog.sciencenet.cn/blog-3376545-1095933.html （有关R的sp包中类型）

假设检验：

https://lrita.github.io/2019/01/28/statistical-hypotheses/

https://lrita.github.io/images/posts/math/%E7%AC%AC%E5%85%AB%E7%AB%A0_%E5%81%87%E8%AE%BE%E6%A3%80%E9%AA%8C.pdf

https://psystat.github.io/chapter/statistics_8.html

https://blog.csdn.net/young_gy/article/details/47865281 （非参数检验）

https://zhuanlan.zhihu.com/p/22987015 （非参数检验）

https://blog.csdn.net/u011517132/article/details/105565678 （关于Wald检验、ML检验、似然比检验）

https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faqhow-are-the-likelihood-ratio-wald-and-lagrange-multiplier-score-tests-different-andor-similar/ （关于Wald检验、ML检验、似然比检验）

https://bbs.pinggu.org/thread-328315-1-1.html （有关Wald检验、ML检验、似然比检验）

https://qastack.cn/stats/6505/likelihood-ratio-test-in-r （似然比检验）

https://www.rdocumentation.org/packages/aod/versions/1.3.1/topics/wald.test （wald检验）

https://www.statisticshowto.com/wp-content/uploads/2016/09/2101f12WaldWithR.pdf （wald检验）

https://rdrr.io/cran/randtests/src/R/runs.test.R （Wald-Wolfowitz游程检验源码）

http://blog.sciencenet.cn/blog-498076-644753.html （分析分类和定序数据的非参数检验）

https://www.datanovia.com/en/lessons/sign-test-in-r/ （SIGN-test）

https://www.datanovia.com/en/lessons/wilcoxon-test-in-r/ （Wilcoxon-test）

https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Cochrans_Q_Test.pdf （Cochran’s Q Test）

相关性分析：

https://www.cnblogs.com/lichunl/p/8953614.html （有关Moran’s I）

http://blog.sina.com.cn/s/blog_4b678be40100o464.html（有关Mantel Test）

https://www.cs.cmu.edu/~neill/papers/jsm2014.pdf （有关Space-Time Interaction）

http://blog.sina.com.cn/s/blog_4b678be40100o464.html（有关Space-Time Interaction）

https://github.com/ricket-sjtu/bi028/wiki/%E5%85%B8%E5%9E%8B%E7%9B%B8%E5%85%B3%E5%88%86%E6%9E%90%EF%BC%88canonical-correlation-analysis,-CCA%EF%BC%89 （典型相关分析）

https://www.biaodianfu.com/acf-pacf.html （时间序列-自相关与偏自相关）

方差检验：

http://biotrainee.com/jmzeng/markdown/ANOVA.html （ANOVA使用）

http://blog.sciencenet.cn/blog-3406804-1190969.html （ANOVA使用）

https://www.datanovia.com/en/lessons/anova-in-r/ （ANOVA使用）

https://www.jianshu.com/p/aa80b6f65399 （ANOVA使用）

https://wap.sciencenet.cn/blog-3406804-1190972.html?mobile=1 （ANOVA非参数检验替代方法）

http://blog.sciencenet.cn/blog-651374-1003788.html （数据顺序对”方差分析ANOVA”的影响）

https://www.datanovia.com/en/lessons/friedman-test-in-r/ （Friedman Test）

https://rdrr.io/cran/rcompanion/man/scheirerRayHare.html （Scheirer Ray Hare test）

回归分析：

https://zhuanlan.zhihu.com/p/52752248 （回归分析数学笔记）

https://zhuanlan.zhihu.com/p/20700731 （有关残差异方差性）

https://wiki.mbalib.com/wiki/%E5%BC%82%E6%96%B9%E5%B7%AE%E6%80%A7（有关残差异方差性）

http://blog.sina.com.cn/s/blog_d81d08a50101gzwy.html （R实现变量优选）

https://flystar233.github.io/2018/09/29/R-regression-analysis/ （R实现变量优选）

https://www.cnblogs.com/pingzeng/p/5073159.html （有关线性回归模型违背基本假定）

https://www.jianshu.com/p/4a44cddd488f （检验、处理异方差性）

https://zhuanlan.zhihu.com/p/27461357 （通过变量筛选解决多重共线性问题）

http://feer.ecnu.edu.cn/_upload/article/files/df/d7/9fe890fc421aa490a209b9e29d03/5ea714f5-b7ac-4c83-a07d-7329847d42b4.pdf （有关内生性）

https://zhuanlan.zhihu.com/p/35518844（有关内生性检验）

http://users.stat.umn.edu/~sandy/courses/8053/handouts/robust.pdf （有关Robust regression）

https://www.jianshu.com/p/4d6f77588b9a （有关Quantile regression）

https://zhuanlan.zhihu.com/p/69816364 （有关含定性变量的回归模型）

https://www.bbsmax.com/A/MAzAoxy59p （分位数回归）

https://cran.r-project.org/web/packages/quantreg/quantreg.pdf （分位数回归）

https://github.com/cloudly/Play-Econometrics-with-R/blob/master/02-cross-section-analysis.md （关于BP检验与white检验）

https://www.zhihu.com/question/53242918 （泊松回归模型和对数线性模型区别）

https://ourcodingclub.github.io/tutorials/mixed-models/ （线性混合效应模型）

https://alexsingleton.files.wordpress.com/2014/09/47-spatial-aurocorrelation.pdf （空间自回归理论基础）

http://stats.lse.ac.uk/lam/SpWgtLag2.pdf (Spatial Lag Model with Time-lagged Effects and Spatial WeightMatrix Estimation理论基础)

聚类分析：

http://blog.sciencenet.cn/blog-3406804-1175517.html

https://www.jianshu.com/p/50cb85285af0

https://www.jianshu.com/p/4f032dccdcef （kmeans聚类算法）

https://rspatial.org/rs/rs.pdf （遥感聚类）

https://www.r-bloggers.com/lang/chinese/619

https://www.cnblogs.com/think90/p/7133753.html

https://blog.csdn.net/sinat_26917383/article/details/51611519

https://blog.csdn.net/brunowuestelle/article/details/46698241 （制作3D可视化图）

https://zhuanlan.zhihu.com/p/23776251 （主成分分析-原理、碎石图）

https://zhuanlan.zhihu.com/p/23795029 （主成分分析-分析、得分、回归）

https://www.jianshu.com/p/2f71bc493042 （主成分分析）

https://zhuanlan.zhihu.com/p/31291210 （因子分析）

http://blog.sciencenet.cn/blog-3406804-1208937.html （典范对应分析CCA）

https://www.jianshu.com/p/bd3820e5a75a （典范对应分析CCA）

http://blog.sciencenet.cn/blog-3406804-1206643.html （CA与CDA分析）

其他：

http://xuewei4d.github.io/

https://www.loyhome.com/

http://bindog.github.io/

http://www2.coe.pku.edu.cn/tpic/file/20170516/20170516145275577557.pdf

https://www3.nd.edu/~mhaenggi/ee87021/Dixon-K-Function.pdf （Ripley’s K）

ps：如果要实现jupyter notebook和R的联动，在安装jupyter notebook基础上（不管是不是在conda里面），R的命令行里面直接使用如下两条命令即可——

1
2
install.packages('IRkernel')
IRkernel::installspec()  # to register the kernel in the current R installation