R语言 为我的实验分析寻找一个适当的工作事后测试,有什么建议吗?[关闭]

tvmytwxo  于 12个月前  发布在  其他
关注(0)|答案(1)|浏览(146)

**已关闭。**此问题不符合Stack Overflow guidelines。目前不接受答案。

这个问题似乎不是关于help center定义的范围内的编程。
2天前关闭。
Improve this question
我之前发布了我的问题,但我编辑了我的问题后没有得到答案,所以我再次发布。
由于正态性问题,我不得不从lmer改为glmmTMB。(emmeans -简单的对比)不知何故不起作用。还有什么建议我可以使用?评论者建议特别是emmeans包,不幸的是,Tukey比较没有回答我的问题,所以我决定使用对比。这将是伟大的,但不与我目前的模型工作。还有什么可能是好的,在我的情况下?
我想检查每个监测年(TIME.SINCE变量)内处理区和对照区(TREATED变量)之间是否存在显著差异。
在这里你可以找到完整的脚本和数据集:https://www.dropbox.com/scl/fo/a54dku8dxwre5wmq8sx2k/h?rlkey=s6rw85iglxl68r1jfea6y2zpg&dl=0

  • 事后脚本:*
library(emmeans)
#time factor!  
#post hoc test is not working with continous variables!!!!!!!! !!!!!! !!!!!
dt2$time.since <- as.factor(dt2$TIME.SINCE)
model2_4_ <- glmmTMB(sqrt(PW_c)~TREATED*time.since+(1|SITE/CODE),data=dt2)
summary(model2_4_)
AIC(model2_4_)
Anova (model2_4_)
#Response: sqrt(PW)
#Chisq Df Pr(>Chisq)
#TREATED 0.9165 1 0.3384012
#time.since 48.6197 11 1.107e-06 ***
#TREATED:time.since 33.8686 10 0.0001944 ***
#simple emmeans contrast
exp2_PW.emm <- emmeans(model2_4_, ~ TREATED * time.since, data=dt2)
contrast(exp2_PW.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")`

字符串

  • 对比结果:*
> contrast(exp2_PW.emm, "consec", simple = "each", combine = TRUE, adjust = "mvt")
 time.since TREATED contrast                    estimate    SE  df t.ratio p.value
 1          .       TREATED1 - TREATED0           nonEst    NA  NA      NA      NA
 2          .       TREATED1 - TREATED0          -0.7130 0.415 254  -1.717      NA
 3          .       TREATED1 - TREATED0          -0.2078 0.323 254  -0.643      NA
 4          .       TREATED1 - TREATED0          -0.4262 0.254 254  -1.680      NA
 5          .       TREATED1 - TREATED0          -0.8181 0.218 254  -3.751      NA
 6          .       TREATED1 - TREATED0          -0.7933 0.211 254  -3.751      NA
 7          .       TREATED1 - TREATED0          -0.4185 0.207 254  -2.022      NA
 8          .       TREATED1 - TREATED0          -0.3727 0.190 254  -1.957      NA
 9          .       TREATED1 - TREATED0          -0.1256 0.164 254  -0.766      NA
 11         .       TREATED1 - TREATED0           0.0197 0.157 254   0.126      NA
 13         .       TREATED1 - TREATED0          -0.0538 0.135 254  -0.398      NA
 21         .       TREATED1 - TREATED0           0.1862 0.131 254   1.419      NA
 .          0       time.since2 - time.since1     0.4418 0.388 254   1.139      NA
 .          0       time.since3 - time.since2     0.3699 0.341 254   1.085      NA
 .          0       time.since4 - time.since3     0.1594 0.246 254   0.647      NA
 .          0       time.since5 - time.since4     0.2679 0.207 254   1.292      NA
 .          0       time.since6 - time.since5     0.1431 0.207 254   0.690      NA
 .          0       time.since7 - time.since6    -0.5211 0.202 254  -2.582      NA
 .          0       time.since8 - time.since7     0.0807 0.191 254   0.423      NA
 .          0       time.since9 - time.since8     0.0373 0.170 254   0.220      NA
 .          0       time.since11 - time.since9   -0.4131 0.143 254  -2.883      NA
 .          0       time.since13 - time.since11   0.4008 0.130 254   3.093      NA
 .          0       time.since21 - time.since13  -0.3651 0.133 254  -2.754      NA
 .          1       time.since2 - time.since1     nonEst    NA  NA      NA      NA
 .          1       time.since3 - time.since2     0.8752 0.358 254   2.447      NA
 .          1       time.since4 - time.since3    -0.0590 0.287 254  -0.206      NA
 .          1       time.since5 - time.since4    -0.1240 0.218 254  -0.568      NA
 .          1       time.since6 - time.since5     0.1679 0.167 254   1.006      NA
 .          1       time.since7 - time.since6    -0.1463 0.165 254  -0.889      NA
 .          1       time.since8 - time.since7     0.1266 0.157 254   0.808      NA
 .          1       time.since9 - time.since8     0.2844 0.135 254   2.099      NA
 .          1       time.since11 - time.since9   -0.2678 0.127 254  -2.106      NA
 .          1       time.since13 - time.since11   0.3272 0.118 254   2.781      NA
 .          1       time.since21 - time.since13  -0.1250 0.109 254  -1.145      NA

Note: contrasts are still on the sqrt scale 
P value adjustment: mvt method for 34 tests

  • 有问题的glmmTMB模型摘要:*
> summary(modEXP2_4<-glmmTMB(sqrt(PW_c)~TREATED*TIME.SINCE+(1|SITE/CODE),data=dt2)) 
 Family: gaussian  ( identity )
Formula:          sqrt(PW_c) ~ TREATED * TIME.SINCE + (1 | SITE/CODE)
Data: dt2

     AIC      BIC   logLik deviance df.resid 
   418.3    443.7   -202.1    404.3      273 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 CODE:SITE (Intercept) 0.06559  0.2561  
 SITE      (Intercept) 0.02366  0.1538  
 Residual              0.20246  0.4500  
Number of obs: 280, groups:  CODE:SITE, 76; SITE, 3

Dispersion estimate for gaussian family (sigma^2): 0.202 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          1.376753   0.166070   8.290  < 2e-16 ***
TREATED1            -0.567529   0.163699  -3.467 0.000526 ***
TIME.SINCE          -0.018044   0.008721  -2.069 0.038551 *  
TREATED1:TIME.SINCE  0.037464   0.010879   3.444 0.000574 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning messages:
1: In sqrt(PW_c) : NaNs produced
2: In sqrt(PW_c) : NaNs produced


谢谢你的帮忙。
我尝试了其他事后测试(例如Tukey配对测试),但它将所有内容与批量中的所有内容进行比较,这不适合我。

wfveoks0

wfveoks01#

你将结果居中,然后取平方根。然而,居中PW会产生许多负值,取sqrt()会产生许多NA值。这可能是你这里问题的原因:

d <- read.csv(
  "c:/Users/Daniel/Downloads/data_2.csv",
  sep = ";",
  dec = ",",
  na.string = "na"
)
d$PW_c <- as.vector(scale(d$PW))
head(sqrt(d$PW_c))
#> [1] NaN NaN NaN NaN NaN NaN
#> Warning message:
#> In sqrt(d$PW_c) : NaNs produced

字符串
我不知道取sqrt的目的是什么。你的数据是高度右偏的:

取sqrt会使偏度变得更糟:

你的结果的自然尺度/范围是什么?也许你可以考虑贝塔或伽马回归。
关于与数字预测器的对比,这应该是可行的:

library(glmmTMB)
library(ggeffects)
d <- read.csv("c:/Users/Daniel/Downloads/data_2.csv", sep = ";", dec = ",", na.string = "na")

d$PW_c <- as.vector(scale(d$PW))
d$ts <- d$TIME.SINCE

m <- glmmTMB(sqrt(PW_c) ~ TREATED * ts + (1 | SITE / CODE), data = d)
hypothesis_test(m, c("TREATED", "ts"), vcov = vcov(m), test = NULL)
hypothesis_test(m, c("TREATED", "ts"), vcov = vcov(m))

相关问题