1.混合模型是否适合您的需求？

str(data)
## 'data.frame':    84 obs. of  6 variables:
##  $Test.ID : int 1 2 3 4 5 6 7 8 9 10 ... ##$ Observer  : Factor w/ 4 levels "Charles","Michelle",..: 1 4 2 4 1 3 2 2 1 2 ...
##  $Relation : Factor w/ 2 levels "Same","Stranger": 1 1 1 1 1 1 1 1 1 1 ... ##$ Aggression: int  4 1 15 2 1 0 2 0 3 10 ...
##  $Tolerance : int 4 34 14 31 4 13 7 6 13 15 ... ##$ Season    : Factor w/ 2 levels "Early","Late": 1 1 1 1 1 1 1 1 1 1 ...

2.哪种概率分布最适合您的数据？

require(car)
## 正在加载所需的包: car
require(MASS)
# 必须为非零的分布
qqp(Aggression, "norm")

# lnorm 表示对数正态
qqp(Aggression, "lnorm")

# qqp需要估计负二项式，泊松和伽玛分布的参数。您可以使用fitdistr函数生成估算值。保存输出并提取每个参数的估计值，如下所示。
fitdistr(rAggression, "Negative Binomial")

qqp(Aggressio, "pois", estimate)

fitdistr(Aggression.t, "gamma")

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

01

02

03

04

3.如何将混合模型拟合到您的数据

3a.如果您的数据是正态分布的

str(starlings)
## 'data.frame':    28 obs. of  5 variables:
##  $Individual : Factor w/ 28 levels "B-40917","B-41205",..: 4 5 6 15 3 16 8 13 20 14 ... ##$ Sex        : Factor w/ 2 levels "F","M": 2 2 2 2 2 1 1 1 1 2 ...
##  $Group : Factor w/ 5 levels "DRT1","MRC1",..: 2 5 5 4 4 4 4 4 4 4 ... ##$ Social.Rank: Factor w/ 2 levels "Breeder","Helper": 2 1 1 1 2 2 2 2 1 2 ...
##  $Mean.Pitch : num 2911 2978 3313 3268 3312 ... summary(lmm) ## Linear mixed model fit by maximum likelihood $'lmerMod'$ ## Formula: Mean.Pitch ~ Sex + Social.Rank + (1 | Group) ## Data: starlings ## ## AIC BIC logLik deviance df.resid ## 389.3 396.0 -189.7 379.3 23 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -2.0559 -0.6272 0.0402 0.5801 2.0110 ## ## Random effects: ## Groups Name Variance Std.Dev. ## Group (Intercept) 0 0 ## Residual 44804 212 ## Number of obs: 28, groups: Group, 5 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 3099.0 82.2 37.7 ## SexM 51.7 81.3 0.6 ## Social.RankHelper -45.0 82.4 -0.5 ## ## Correlation of Fixed Effects: ## (Intr) SexM ## SexM -0.630 ## Scl.RnkHlpr -0.668 0.106 让我们看看结果。首先，我们获得一些模型拟合的度量，包括AIC，BIC，对数似然度和偏差。然后我们得到由随机效应解释的方差估计。这个数字很重要，因为如果它与零没有区别，那么您的随机效应可能并不重要，您可以继续进行常规的线性模型建立。接下来，我们将对固定效应进行估算，带有标准误差。这些信息可能足以满足您的需求。一些期刊将这些模型的结果报告为带有置信区间的效应大小。当然，当我查看固定效应估算值时，我已经可以看出，性别和社会地位之间的平均音高没有差异。但是有些期刊希望您报告p值。 如果您想要一些p值，则需要使用Anova函数。 ## Analysis of Deviance Table (Type II Wald chisquare tests) ## ## Response: Mean.Pitch ## Chisq Df Pr(>Chisq) ## Sex 0.4 1 0.52 ## Social.Rank 0.3 1 0.58 Anova函数进行了Wald检验，该检验告诉我们我们对性别和社会地位对音高的影响的估计p值。 拟合线性混合模型时，可能会遇到一种复杂情况。R可能会有“无法收敛”错误，通常将其表述为“没有收敛就达到了迭代限制”。这意味着您的模型有太多因素，样本量不够大，无法拟合。然后，您应该做的是从模型中删除固定效果和随机效果，然后进行比较以找出最合适的效果。一次删除固定效果和随机效果。保持固定效果不变，并一次删除一个随机效果，然后找出最合适的效果。然后保持随机效果不变，并一次删除固定效果。在这里，我只有一种随机效果， anova(noranklmm, nosexlmm, nofixedlmm) ## Data: starlings ## Models: ## nofixedlmm: Mean.Pitch ~ 1 + (1 | Group) ## noranklmm: Mean.Pitch ~ Sex + (1 | Group) ## nosexlmm: Mean.Pitch ~ Social.Rank + (1 | Group) ## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) ## nofixedlmm 3 386 390 -190 380 ## noranklmm 4 388 393 -190 380 0.48 1 0.49 ## nosexlmm 4 388 393 -190 380 0.00 0 1.00 请注意，该方差分析函数与我们用来评估模型中固定效果的重要性的方差分析函数不同。方差分析函数用于比较模型。p值表明模型之间没有明显的重要差异。我们还可以比较AIC值，请注意，具有最低AIC值的模型是完全没有固定影响的模型，这符合我们的理解，即性别和社会地位对歌曲的音调没有影响。无论采用哪种方法，请务必在稿件中报告用于选择最佳模型的p值或AIC值。 3b.如果您的数据不是正态分布的 您会看到，用于估计模型中影响大小的REML和最大似然法做出了不适用于数据的正态假设，因此您必须使用其他方法进行参数估计。问题在于，存在许多替代的估算方法，每种估算方法都使用不同的R包运行，并且很难确定哪种方法合适。 首先，我们需要测试是否可以使用惩罚拟似然（PQL）。PQL是一种灵活的技术，可以处理非正常数据，不平衡设计和交叉随机效应。但是，如果您的因变量符合离散计数分布（例如泊松或二项式）且均值小于5，或者您的因变量为二元变量，则会产生偏差估计。 Aggression变量适合对数正态分布，该分布不是离散分布。这意味着我们可以继续使用PQL方法。但是在继续之前，让我们回到转变为正态的问题。 将分布设置为对数正态，我们将族设置为高斯，并将链接设置为log。 ## lmList summary(PQL) ## Linear mixed-effects model fit by maximum likelihood ## Data: recog ## AIC BIC logLik ## NA NA NA ## ## Random effects: ## Formula: ~1 | Observer ## (Intercept) ## StdDev: 0.3312 ## ## Formula: ~1 | Test.ID %in% Observer ## (Intercept) Residual ## StdDev: 0.5295 7.128 ## ## Variance function: ## Structure: fixed weights ## Formula: ~invwt ## Fixed effects: Aggression.t ~ Relation + Season ## Value Std.Error DF t-value p-value ## (Intercept) 1.033 0.5233 55 1.974 0.0535 ## RelationStranger 1.210 0.4674 55 2.589 0.0123 ## SeasonLate -1.333 0.5983 23 -2.228 0.0359 ## Correlation: ## (Intr) RltnSt ## RelationStranger -0.855 ## SeasonLate -0.123 0.000 ## ## Standardized Within-Group Residuals: ## Min Q1 Med Q3 Max ## -4.86916 -0.29958 -0.08012 0.14280 5.93336 ## ## Number of Observations: 84 ## Number of Groups: ## Observer Test.ID %in% Observer ## 4 28 该模型表明季节对侵略性有影响，也就是说，在菌落周期后期收集的黄蜂比早期收集的黄蜂侵略性小。这也表明黄蜂之间的关系有影响。他们相对于巢友更可能对陌生人有侵略性。我将这些统计数据与估计值，标准误，t值和p值一起报告。 那么，如果您的因变量的平均值小于5，或者您有一个二元因变量，而您不能使用PQL，该怎么办？在这里，您可以使用两种选择：拉普拉斯（Laplace）近似 和马尔可夫链蒙特卡罗算法（MCMC）。拉普拉斯近似值最多可以处理3个随机效果。除此之外，您还必须使用MCMC。 让我们从一个可以使用拉普拉斯逼近的例子开始。我们将使用学生在学校的学习情况的数据。出于本示例的目的，我将仅将数据子集化为几个感兴趣的变量，并将“ repeatgr”变量简化为二元因变量。 str(bdf) ## 'data.frame': 2287 obs. of 4 variables: ##$ schoolNR: Factor w/ 131 levels "1","2","10","12",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $Minority: Factor w/ 2 levels "N","Y": 1 2 1 1 1 2 2 1 2 2 ... ##$ ses     : num  23 10 15 23 10 10 23 10 13 15 ...
##  $repeatgr: Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 2 1 ... 假设我们要找出是否属于少数民族和社会经济地位会影响学生复读成绩的可能性。我们的因变量是“ repeatgr”，指示学生是否重复了成绩。少数族身份是二元Y / N类别，社会经济地位由“ ses”表示，其数字范围为10至50，其中50是最富有的。我们的随机因素是“ schoolNR”，它代表从中采样学生的学校。因为因变量是二元的，所以我们需要具有二项式分布的广义线性混合模型，并且由于我们的随机效应少于五个，因此可以使用Laplace近似 。 严格来说，Laplace近似 是一种称为Gauss-Hermite正交（GHQ）的参数估算方法的特例，需要进行一次迭代。GHQ由于重复迭代而比Laplace更为准确，但在第一次迭代后变得不那么灵活，因此您只能将它用于一个随机效果。我们唯一的随机效应是“ schoolNR”。 summary(GHQ) ## Generalized linear mixed model fit by maximum likelihood (Adaptive ## Gauss-Hermite Quadrature, nAGQ = 25) $glmerMod$ ## Family: binomial ( logit ) ## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR) ## Data: bdf ## ## AIC BIC logLik deviance df.resid ## 1672.8 1701.4 -831.4 1662.8 2282 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -0.964 -0.407 -0.314 -0.221 5.962 ## ## Random effects: ## Groups Name Variance Std.Dev. ## schoolNR (Intercept) 0.264 0.514 ## Number of obs: 2287, groups: schoolNR, 131 ## ## Fixed effects: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.45194 0.20763 -2.18 0.03 * ## MinorityY 0.47957 0.47345 1.01 0.31 ## ses -0.06205 0.00798 -7.78 7.5e-15 *** ## MinorityY:ses 0.01196 0.02289 0.52 0.60 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Correlation of Fixed Effects: ## (Intr) MnrtyY ses ## MinorityY -0.400 ## ses -0.905 0.369 ## MinortyY:ss 0.299 -0.866 -0.321 Laplace <- glmer(repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR), data = bdf, family = binomial(link = "logit")) # Contrast to the Laplace approximation ## Warning: Model failed to converge with max|grad| = 0.0458484 (tol = 0.001) summary(Laplace) ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) $glmerMod$ ## Family: binomial ( logit ) ## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR) ## Data: bdf ## ## AIC BIC logLik deviance df.resid ## 1672.8 1701.5 -831.4 1662.8 2282 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -0.960 -0.407 -0.316 -0.222 5.950 ## ## Random effects: ## Groups Name Variance Std.Dev. ## schoolNR (Intercept) 0.258 0.508 ## Number of obs: 2287, groups: schoolNR, 131 ## ## Fixed effects: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.45456 0.20611 -2.21 0.027 * ## MinorityY 0.48005 0.47121 1.02 0.308 ## ses -0.06191 0.00791 -7.83 4.9e-15 *** ## MinorityY:ses 0.01194 0.02274 0.53 0.600 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Correlation of Fixed Effects: ## (Intr) MnrtyY ses ## MinorityY -0.400 ## ses -0.906 0.369 ## MinortyY:ss 0.299 -0.866 -0.321 您可以看到在这种情况下，Laplace和GHQ之间没有重要区别。两者都表明，社会经济课对学生复读成绩的可能性有非常显着的影响，尽管即使采用logit变换，我们也可以看到影响量很小。但是，使用此方法时，还需要考虑其他一些因素。当用于过度分散的数据时，即合并的残差远大于残差自由度时，它变得不准确。使用此方法时，应检查模型以确保数据不会过度分散。 ## n×n方差-协方差矩阵中方差参数的数量 vpars <- function(m) { nrow(m) * (nrow(m) + 1)/2 } # 接下来计算剩余自由度 rdf <- nrow(model.frame(model)) - model.df # 提取皮尔逊残差 prat <- Pearson.chisq/rdf # 生成一个p值。如果小于0.05，则数据过于分散。 这些学校数据并不分散。如果是，可以将过度分散建模为随机效应，每个观察值具有一个随机效应水平。在这种情况下，我将使用学生ID作为随机效果变量。 我们继续讨论泊松数据的均值太小或因变量是分类的情况，并且我们有五个或五个以上随机效应。考虑有关大麦农民的这些数据。想象一下，要使大麦收成产生利润，大麦收成的收入必须大于140。 farmers <- farmers$c(1:5, 8)$ str(farmers) ## 'data.frame': 102 obs. of 6 variables: ##$ year    : int  1901 1901 1901 1901 1902 1902 1902 1902 1902 1902 ...
##  $farmer : Factor w/ 18 levels "Allardyce","Dooley",..: 10 5 3 18 10 5 18 17 4 12 ... ##$ place   : Factor w/ 16 levels "Arnestown","Bagnalstown",..: 3 16 14 11 3 16 11 4 8 6 ...
##  $district: Factor w/ 4 levels "CentralPlain",..: 2 2 1 1 2 2 1 1 4 4 ... ##$ gen     : Factor w/ 2 levels "Archer","Goldthorpe": 1 1 1 1 1 1 1 1 1 1 ...
##  \$ profit  : num  1 1 1 1 1 1 1 1 1 1 ...

##  Iterations = 3001:12991
##  Thinning interval  = 10
##  Sample size  = 1000
##
##  DIC: 76.61
##
##  G-structure:  ~year
##
##      post.mean l-95% CI u-95% CI eff.samp
## year      18.2    0.266     68.6     18.1
##
##                ~farmer
##
##        post.mean l-95% CI u-95% CI eff.samp
## farmer      1.93   0.0789     6.68     99.7
##
##                ~place
##
##       post.mean l-95% CI u-95% CI eff.samp
## place      1.54   0.0711     4.84      173
##
##                ~gen
##
##     post.mean l-95% CI u-95% CI eff.samp
## gen      12.1   0.0874     28.6     1000
##
##                ~district
##
##          post.mean l-95% CI u-95% CI eff.samp
## district      1.78   0.0639     5.63     1000
##
##  R-structure:  ~units
##
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
##
##  Location effects: profit ~ 1
##
##             post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept)      3.47    -1.16    10.75      220  0.14

结束：了解您的数据

ggplot(recog, aes(x = Aggression)+ geom_density() +

Trace(MCMC, log = TRUE)

Trace(MCMC2, log = TRUE)

# 将两个模型的估计值和置信区间放在一起
rbind (covariances, Gcovariances)
# 创建一个数据框架，其中包含模型和随机效应的因素
data.frame(coint, model = factor(el1", "model2"), eac),
random.effect = dimnames(i)
## Warning: some row.names duplicated: 6,7,8,9,10 --> row.names NOT used
# 将估计和置信区间绘制为按模型分组的箱形图
ggplot(conf.int+ geom_crossbar(aes(y.95..CI,
y.95..CI= model= "dodge")

|
1月前
【R语言实战】——带有高斯新息的金融时序的GARCH模型拟合预测及VAR/ES风险度量
【R语言实战】——带有高斯新息的金融时序的GARCH模型拟合预测及VAR/ES风险度量
41 2
|
1月前
【R语言实战】——带有新息为标准学生t分布的金融时序的GARCH模型拟合预测
【R语言实战】——带有新息为标准学生t分布的金融时序的GARCH模型拟合预测
33 1
|
1月前
【R语言实战】——Logistic回归模型
【R语言实战】——Logistic回归模型
38 1
|
1月前
|

R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码2
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码
40 2
|
1月前
|

R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
60 1
|
1月前
|

【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
47 2
|
1月前
|

R语言逻辑回归logistic模型ROC曲线可视化分析2例：麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例：麻醉剂用量影响、汽车购买行为
53 2
|
1月前
|

305 3
|
1月前

37 3
|
1月前
|
Web App开发 数据可视化 数据挖掘

151 2