我正在用R编写循环或函数,但我仍然没有真正理解如何做到这一点。目前,我需要编写一个循环/函数(不确定哪个更好),以便在排列测试或随机化中创建混合模型公式的几个结果
示例数据集如下所示:
dataset <- read.table(text =
"ID A_2 B_2 C_2 A_1 B_1 C_1 chkgp
M1 10 20 60 30 54 33 Treatment
M1 20 50 40 33 31 44 Placebo
M2 40 80 40 23 15 66 Placebo
M2 30 90 40 67 67 66 Treatment
M3 30 10 20 22 89 77 Treatment
M3 40 50 30 44 50 88 Placebo
M4 40 30 40 42 34 99 Treatment
M4 30 40 50 33 60 80 Placebo",header = TRUE, stringsAsFactors = FALSE)
每次运行模型时,我都会随机处理chkgp变量,并使用以下代码
mod1<-summary(lmerTest::lmer(A_2~B_1+sample(chkgp)+(1|ID),data = dataset))
mod1
P_value= 2 * (1 - pnorm(abs(mod1$coefficients[3, 4])))
P_value
混洗1次后的结果模型1
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: A_2 ~ B_1 + sample(chkgp) + (1 | ID)
Data: dataset
REML criterion at convergence: 44.7
Scaled residuals:
Min 1Q Median 3Q Max
-0.7792 -0.4441 0.1185 0.3893 0.7734
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 122.52 11.069
Residual 10.15 3.186
Number of obs: 8, groups: ID, 4
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 42.87450 6.49941 4.06743 6.597 0.00258 **
B_1 -0.27492 0.09149 7.97033 -3.005 0.01702 *
sample(chkgp)Treatment 1.74313 4.71095 8.33366 0.370 0.72060
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) B_1
B_1 -0.433
smpl(chkg)T 0.164 -0.748
[1] 0.7113693
问题1:我需要找出比较和检查实际P值是否相同的方法,即使在混洗chkgp变量1000之后也是如此
问题2:我需要编写循环运行模型1000次,每次都需要对chkgp变量进行混洗。
发布于 2019-05-30 01:52:15
dataset <- read.table(text =
"ID A_2 B_2 C_2 A_1 B_1 C_1 chkgp
M1 10 20 60 30 54 33 Treatment
M1 20 50 40 33 31 44 Placebo
M2 40 80 40 23 15 66 Placebo
M2 30 90 40 67 67 66 Treatment
M3 30 10 20 22 89 77 Treatment
M3 40 50 30 44 50 88 Placebo
M4 40 30 40 42 34 99 Treatment
M4 30 40 50 33 60 80 Placebo",header = TRUE, stringsAsFactors = FALSE)
storeP <- list(rep(NA, 1000))
for(i in 1:1000) {
#dataset <- transform(dataset, chkgp = sample(chkgp))
dataset$chkgp <- ifelse(runif(nrow(dataset)) > 0.5, "Treatment", "Placebo")
mod1<-summary(lmerTest::lmer(A_2~B_1+sample(chkgp)+(1|ID),data = dataset))
P_value= 2 * (1 - pnorm(abs(mod1$coefficients[3, 4])))
storeP[[i]] <- P_value
}
var(unlist(storeP)) == 0
我有两个“数据混洗”的选项:
chkgp:transform(dataset, chkgp = sample(chkgp))
列的
以0.5阈值随机分配“治疗”和“安慰剂”给变量:ifelse(runif(nrow(dataset)) > 0.5, "Treatment", "Placebo"
:
https://stackoverflow.com/questions/56365728
复制相似问题