我想执行一个ks.test来比较两个发行版。我在考虑使用Kolmogorov-Smirnov检验,但问题是这两个分布都必须加权。你知道怎么做吗?
以下是我的数据:
library(tidyverse)
my_data_2018 <- tibble(Var = c(900, 1500, 350, 1200, 750, 100,125,250),
my_weights_2018 = c(2.2, 3.1, 8.2, 4.2, 5.3, 6.8, 12, 25))
my_data_2019 <- tibble(Var = c(32, 21, 21, 900, 1500, 350, 1200, 750, 100,125,250,300),
my_weights_2019 = c(2.2, 3.1, 8.2, 2.2, 3.1, 8.2, 4.2, 5.3, 6.8, 12, 25, 1))
我找到这段代码来创建一个新的ks_weighted函数,但我不确定它如何处理我的示例数据
ks_weighted <- function(vector_1,vector_2,weights_1,weights_2){
F_vec_1 <- ewcdf(vector_1, weights = weights_1, normalise=FALSE)
F_vec_2 <- ewcdf(vector_2, weights = weights_2, normalise=FALSE)
xw <- c(vector_1,vector_2)
d <- max(abs(F_vec_1(xw) - F_vec_2(xw)))
## P-VALUE with NORMAL SAMPLE
# n_vector_1 <- length(vector_1)
# n_vector_2<- length(vector_2)
# n <- n_vector_1 * n_vector_2/(n_vector_1 + n_vector_2)
# P-VALUE EFFECTIVE SAMPLE SIZE as suggested by Monahan
n_vector_1 <- sum(weights_1)^2/sum(weights_1^2)
n_vector_2 <- sum(weights_2)^2/sum(weights_2^2)
n <- n_vector_1 * n_vector_2/(n_vector_1 + n_vector_2)
pkstwo <- function(x, tol = 1e-06) {
if (is.numeric(x))
x <- as.double(x)
else stop("argument 'x' must be numeric")
p <- rep(0, length(x))
p[is.na(x)] <- NA
IND <- which(!is.na(x) & (x > 0))
if (length(IND))
p[IND] <- .Call(stats:::C_pKS2, p = x[IND], tol)
p
}
pval <- 1 - pkstwo(sqrt(n) * d)
out <- c(KS_Stat=d, P_value=pval)
return(out)
}
发布于 2020-06-22 22:28:15
您可以使用FastDR中的ks
函数来获得加权的Kolmogorov-Smirnov统计量,然后计算双边检验的p值。
FastDR是由UPenn的教授Greg Ridgeway撰写的。我认为这个方法是可靠的。
首先,从GitHub获取包。
if(!require(devtools)) install.packages("devtools")
library(devtools)
install_github("gregridgeway/survey","patch-1")
install_github("gbm-developers/gbm3")
install_github("gregridgeway/fastDR")
要根据您的数据计算KS测试统计数据,您可以使用以下代码:
library(fastDR)
Var <- c(my_data_2018$Var,my_data_2019$Var)
weights <- c(my_data_2018$my_weights_2018,my_data_2019$my_weights_2019)
treat <- c(rep(0,nrow(my_data_2018)),rep(1,nrow(my_data_2019)))
ks.test.stats <- ks(x=Var,z=treat,w=weights)
您将获得以下KS测试统计信息:
## KS Test Statistics
ks.test.stats
[1] 0.1660517
最后,要计算双边测试的p值,可以使用以下代码:
### "Two-Sided" P-Value for KS Tests Statistics from R
p.value <- 1-.Call(stats:::C_pSmirnov2x,
STAT=ks.test.stats,
length(my_data_2018$Var),
length(my_data_2019$Var))
p.value
> ks.test.stats
[1] 0.1660517
注意:如果您在下载软件包时遇到问题,请尝试更改R版本。FastDR的依赖项在R版本4.00中似乎不能很好地工作。我使用R版本3.63来执行此计算。
https://stackoverflow.com/questions/62513941
复制相似问题