x<-c(221.60112, 305.217725, 295.251684, 371.3313, 397.452722, 437.212724)
group<-factor(c(rep(1,3), rep(2, 3)))
df<-data.frame(x, group)
library(coin)
oneway_test(x~group, data=df, distribution="exact")
wilcox_test(x~group, data=df, distribution="exact", conf.int=T)

# Обратите внимание, что границы 95% доверительного интервала для медианной разности не определены.

#Понизим доверительный уровень:
wilcox_test(x~group, data=df, distribution="exact", conf.int=T, conf.level = 0.9)
# Нижняя границы 90% ДИ определена, верхняя - бесконечность.

# А теперь фокус-покус:
df2<-df
df2$x[4:6]<-df$x[4:6]*100000

#Этим мы увеличили значения всех вариант выборки ?2 в сто тысяч раз, то есть смоделировлан случай, когда и безо всякого статанализа колоссальные различия в параметрах положения распределений очевидны.

# Парадоксально, но уровень значимости точных тестов не поменялся:

oneway_test(x~group, data=df2, distribution="exact")
wilcox_test(x~group, data=df2, distribution="exact")

# А все потому, что тестовая статистика при перестановке может принимать всего лишь 20 значений:
choose(6, 3)

# Так что использование рандомизации для решения данной конкретной задачи - так себе идея.     
# А вот бутсреп работает:

library(confintr)
ci_mean_diff(x=df$x[1:3], y=df$x[4:6], type="boot", boot_type="bca", R=99999)

# Более того, даже 99,9% ДИ для разности средних не включает 0:
ci_mean_diff(x=df$x[1:3], y=df$x[4:6], type="boot", boot_type="bca", probs =c(0.0005, 0.9995), R=99999)

# То есть гипотеза о равенстве математических ожиданий откланяется на уровне значимости 0,001, чего, если я правильно понял, и желает ТС.

x<-c(221.60112, 305.217725, 295.251684, 371.3313, 397.452722, 437.212724)
group<-factor(c(rep(1,3), rep(2, 3)))
df<-data.frame(x, group)

t.test(x~group, data=df, var.eq=T)

t.test(x~group, data=df)

x<-c(221.60112, 305.217725, 295.251684, 371.3313, 397.452722, 437.212724)
R<-9999
absdiff<-abs(mean(x[1:3])-mean(x[4:6]))
res<-numeric(R)
for (i in 1:R)
{
perm<-sample(x)
permdiff<-abs(mean(perm[1:3])-mean(perm[4:6]))
if(permdiff>=absdiff) res [i]<-1
}
p<-(sum(res)+1)/(R+1)
p

x<-c(221.60112, 305.217725, 295.251684)
y<-c(371.3313, 397.452722, 437.212724)

diff<-274-402
res<-rep(0, 10000)
for (i in 1:10000)
{
x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(y), sd(y))
ci<- t.test(x_sim, y_sim)$conf.int
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
#Доля попаданий:
sum(res)/10000

diff<-274-402
res<-rep(0, 10000)
for (i in 1:10000)
{
x_sim<-c(rnorm(2, mean(x), sd(x)*1/3), rnorm(1, mean(x), sd(x)*2/3*2.5))
y_sim<-c(rnorm(2, mean(y), sd(y)*1/3), rnorm(1, mean(y), sd(y)*2/3*2.5))
ci<- t.test(x_sim, y_sim)$conf.int
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
sum(res)/10000

library(confintr)
diff<-274-402
res<-rep(0, 1000)
for (i in 1:1000)
{
x_sim<-c(rnorm(2, mean(x), sd(x)*1/3), rnorm(1, mean(x), sd(x)*2/3*2.5))
y_sim<-c(rnorm(2, mean(y), sd(y)*1/3), rnorm(1, mean(y), sd(y)*2/3*2.5))
ci<- ci_mean_diff(x_sim, y_sim, type="boot", boot_type="bca")$interval
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
sum(res)/1000

x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(y), sd(y))

## Создадим данные
a <- c(221.60112, 305.217725, 295.251684)
b <- c(371.3313, 397.452722, 437.212724)
ab <- c(a, b)
n <- length(ab)
g <- rep(c('a', 'b'), time = c(length(a), length(b)))

## Трансформация к диапазону (0; 1)
vec <- scales::rescale_max(ab,
from = c(
1/max(ab)/n,
max(ab)+max(ab)/n
)
)

## Бета-регрессия
fit <- betareg::betareg(vec~g)
emm <- emmeans::emmeans(fit, ~g)

## Средние и ДИ для групп
res <- data.frame(emm, row.names = 'g')[c('emmean', 'asymp.LCL', 'asymp.UCL')]
cev <- apply(res, 2, function(x) scales::rescale_max(x,
to = c(
1/max(ab)/n,
max(ab)+max(ab)/n
),
from = c(0, 1)))

## Разница средних и ДИ
ctr <- data.frame(confint(pairs(emm)))[c('estimate', 'asymp.LCL', 'asymp.UCL')]
rtc <- apply(ctr, 2, function(x) scales::rescale_max(x,
to = c(
1/max(ab)/n,
max(ab)+max(ab)/n
),
from = c(0, 1)))

## Расчёт средних в группах
set.seed(32167)
mean.a <- Hmisc::smean.cl.boot(a, B = 9999)
mean.b <- Hmisc::smean.cl.boot(b, B = 9999)
means <- rbind(cev, mean.a, mean.b)
colnames(means) <- c('mu', 'lwr', 'upr')

## Расчёт разницы средних и ДИ
set.seed(32167)
ci.bca <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'bca')
ci.per <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'perc')
ci.nor <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'norm')
ci.bas <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'basic')
fun_extr <- function(x) { #Для единообразного извлечения
estimate <- x[['estimate']]
interval <- x[['interval']]
result <- c(estimate, interval)
return(result)
}
diffs <- rbind(
beta = rtc,
bca = fun_extr(ci.bca),
perc = fun_extr(ci.per),
norm = fun_extr(ci.nor),
basic = fun_extr(ci.bas)
)
colnames(diffs) <- c('mu', 'lwr', 'upr')

## Выводим результат
summary(fit)
print(means, digits = 5)
print(diffs, digits = 5)

> Call:
betareg::betareg(formula = vec ~ g)

Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7892 -1.0236 0.2118 0.9816 1.5627

Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.1480 0.1607 0.921 0.357
gb 1.1632 0.2525 4.606 4.1e-06 ***

Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 50.91 29.16 1.746 0.0808 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type of estimator: ML (maximum likelihood)
Log-likelihood: 8.148 on 3 Df
Pseudo R-squared: 0.7817
Number of iterations: 68 (BFGS) + 2 (Fisher scoring)
> mu lwr upr
a 273.88 233.94 313.82
b 401.80 369.21 434.40
mean.a 274.02 221.60 305.22
mean.b 402.00 371.33 437.21
> mu lwr upr
beta -127.92 -179.45 -76.395
bca -127.98 -191.06 -86.850
perc -127.98 -180.40 -78.143
norm -127.98 -180.30 -74.747
basic -127.98 -177.81 -75.553

##Дополнительные функции #fix
fun_int <- function(x) {
int <- x[['interval']]
return(int)
}
fun_res <- function(data, diff) {
len <- length(data)
res <- vapply(X = 1:n, FUN.VALUE = logical(1),
function(i) {
cond <- diff>data[[i]][1] & diff<data[[i]][2]
return(cond)
}
)
pr <- sum(res)/len
sums <- rowMeans(simplify2array(data))
res <- list(low = sums[1], up = sums[2], l = abs(sums[1]-sums[2]), prob = pr)
return(res)
}

x <- c(221.60112, 305.217725, 295.251684)
y <- c(371.3313, 397.452722, 437.212724)
mu_x <- mean(x)
mu_y <- mean(y)
len_x <- length(x)
sd_x <- sd(x)
sd_y <- sd(y)
len_y <- length(y)
diff <- mu_x-mu_y
len <- len_x+len_y
g <- rep(c('x', 'y'), time = c(len_x, len_y))
idx <- c('asymp.LCL', 'asymp.UCL')

n <- 1000
res <- vector(mode = 'logical', length = n)

set.seed(32167)
x_sim <- replicate(n,
c(
rnorm(round(len_x*2/3), mu_x, sd_x*1/3),
rnorm(round(len_x*1/3), mu_x, sd_x*2/3*2.5)
)
)
y_sim <- replicate(n,
c(
rnorm(round(len_x*2/3), mu_y, sd_y*1/3),
rnorm(round(len_x*1/3), mu_y, sd_y*2/3*2.5)
)
)
xy_sim <- rbind(x_sim, y_sim)
vec_sim <- apply(xy_sim, 2, function(i) scales::rescale_max(i,
from = c(
1/max(i)/len,
max(i)+max(i)/len
)
))

ci_beta <- lapply(
X = 1:n, FUN = function(i) {
fit <- betareg::betareg(vec_sim[,i]~g)
ci <- as.numeric(confint(pairs(emmeans::emmeans(fit, ~g)))[idx])
ic <- scales::rescale_max(ci,
to = c(
1/max(xy_sim[,i])/len,
max(xy_sim[,i])+max(xy_sim[,i])/len
),
from = c(0, 1)
)
return(ic)
})
ci_bca <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'bca', R = n)
ci <- fun_int(fit)
return(ci)
})
ci_perc <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'perc', R = n)
ci <- fun_int(fit)
return(ci)
})
ci_norm <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'norm', R = n) #fix
ci <- fun_int(fit)
return(ci)
})
ci_basic <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'basic', R = n) #fix
ci <- fun_int(fit)
return(ci)
})
vct <- ls(pattern = 'ci_')

res_ci <- t(simplify2array(lapply(1:length(vct), function(i) {
#do.call(fun_res, list(data = get(vct[i]), diff = get('diff')))
fun_res(data = get(vct[i]), get('diff')) #fix
})))
row.names(res_ci) <- vct
print(res_ci, digits = 5) # fix

low up l prob
ci_basic -173.1 -81.769 91.328 0.885
ci_bca -176.03 -81.595 94.435 0.773
ci_beta -172.44 -79.885 92.552 0.866
ci_norm -174 -80.332 93.666 0.876
ci_perc -172.38 -81.104 91.271 0.814

## Сумма
do.call(sum, list(c(1, 2, 4, 1, 2), na.rm = TRUE)) #Правильно
lapply(c(1, 2, 4, 1, 2), sum) #Неправильно

##Структура операций
L <- list(c(1,2,3), c(4,5,6))
L
[[1]]
[1] 1 2 3

[[2]]
[1] 4 5 6

lapply(L, sum)
[[1]]
[1] 6

[[2]]
[1] 15

list( sum(L[[1]]) , sum(L[[2]]))
[[1]]
[1] 6

[[2]]
[1] 15

do.call(sum, L)
[1] 21
sum(L[[1]], L[[2]])
[1] 21

##Передача аргументов в качестве листа
args <- list(n = 10, mean = 50, sd = 10) #Т-баллы
do.call(rnorm, args) #эквивалентно rnorm(args$n, args$mean, args$sd)

replicate(10, rnorm(10)) #В каждом столбце числа разные
do.call(replicate, list(10, rnorm(10))) #В каждом столбце числа одинаковые

##Вызов функции по её названию
do.call('colnames', list(iris))

##Комбинация списков
rbind(L) #do.call(rbind, list(L))
Reduce('rbind',x = L) #do.call(rbind, L)

lapply(c(1, 2, 4, 1, 2), sum) #Неправильно

CWTest<-function(x,y)

{
m<-length(x)
n<-length(y)

mean_X<-mean(x)
mean_Y<-mean(y)

Var_X<-var(x)
Var_Y<-var(y)

Cmn<-sqrt(m*n)*(mean_X-mean_Y)/sqrt(n*Var_X+m*Var_Y)

return(c(TEST=Cmn,p_value=2*pnorm(-abs(Cmn))))

}

> CWTest(x,y)

         TEST       p_value
-3.9267209282  0.0000861118

res<-numeric(10000)
for(i in 1:10000)
{x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(x), sd(y))
res[i]<-CWTest(x_sim, y_sim)[2]
}
hist(res)

Tobs<-CWTest(x, y)[1]
res<-rep(0, 99999)
for(i in 1:99999)
{x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(x), sd(y))
Tsim<-CWTest(x_sim, y_sim)[1]
if(abs(Tsim)>=abs(Tobs))res[i]<-1
}
p=(sum(res)+1)/100000
p

res[i]<-CWTest(x_sim, y_sim)[2]

hist(replicate(500000,CWTest(rnorm(10,mean(x),sd(x)),rnorm(10,mean(x),sd(y)))[1]),col='azure3',main="",xlab="")

hist(replicate(500000,CWTest(rnorm(200),rnorm(200))[2]),col='azure3',main="",xlab="")

R<-100#Максимальный выборочный объем
N<-100000#Количество Монте-Карловых выборок
tab<-rep(0,R-2)

for(i in 3:R)# Для выборочных объемов от 3 до 100 с шагом 1
{
res<-replicate(N,CWTest(rnorm(i), rnorm(i))[1])
tab[i-2]<-length(res[abs(res)>1.96])/N
}

plot(3:R,tab,type="l",main="Ошибка I рода для CWTest()")
abline(h=.05,col='red')

d<-1
p<-replicate(1000, t.test(rnorm(2, 0, 1), rnorm(1, d, 1),var.eq=T)$p.value)
delta<-rep(d, 1000)
res<-data.frame(delta, p)
mean_p1<-mean(p)
mean_p<-rep(NA, 99)
for (i in 1:99)
{
d<-d+1
p<-replicate(1000, t.test(rnorm(2, 0, 1), rnorm(1, d, 1),var.eq=T)$p.value)
mean_p[i]<-mean(p)
delta<-rep(d, 1000)
res_<-data.frame(delta, p)
res<-rbind(res, res_)
}

boxplot(p~delta, data=res, range=0)
abline(0.05,0, lwd=2, col="red")