# 基于R语言混合效应模型（mixed model）案例研究-2

https://developer.aliyun.com/article/1489317

## 3b.如果您的数据不是正态分布的

Aggression变量适合对数正态分布，该分布不是离散分布。这意味着我们可以继续使用PQL方法。但是在继续之前，让我们回到转变为正态的问题。

##     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

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 ...

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

## 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，则数据过于分散。

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语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
47 2
|
1月前
|

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

306 3
|
1月前

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

151 2
|
1月前
|

R语言分类回归分析考研热现象分析与考研意愿价值变现
R语言分类回归分析考研热现象分析与考研意愿价值变现
45 2
|
1月前
|

R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
33 1