Вычисление показателей описательной статистики
Описание
Настраиваемая функция для вычисления показателей описательной статистики и оценки эмпирических распределений. Функция позволяет учитывать несколько факторов в качестве группирующих признаков. Результаты обработки представлены классом “таблица” (data.frame).
Использование
descrstats2 <- function(x, grx, nmin = 3, Mean = TRUE, SE = TRUE, Mean_ng = TRUE, SE_ng = TRUE, Me = TRUE, Min = TRUE, Max = TRUE, Range = TRUE, As = FALSE, Ex = FALSE, Q1 = TRUE, Q3 = TRUE, IQR = TRUE, SD = TRUE, Cv = TRUE, AD = FALSE, SF = FALSE, srclevs = FALSE)
Аргументы
x исходные данные (зависимые признаки с непрерывным характером распределений). grx группирующие признаки. nmin минимальное количество вариант в градации фактора (по умолчанию nmin = 3). … AD критерий Андерсона-Дарлинга. При AD = TRUE в таблицу включается результат соответствующего тестирования и оценка уровня статистической значимости. SF критерий Шапиро-Франциа. При SF = TRUE в таблицу включается результат соответствующего тестирования и оценка уровня статистической значимости.
Исходная функция
#Функция усечения количества вариант в градации фактора n >= ingroup
n_group_restricted <- function(c,ingroup) #c - таблица, где 1 столбец - зависимый признак,
#2-й столбец - группирующий признак;
#ingroup - минимальное количество вариант в градации
{local({
#b - фактор c урезанным кол-вом строк "NA"
#c - фактор без урез.строк "NA", ingroup - количество вариант в градации
names(c) <- c("depvar","factor")
#Урезаем базу с фактором и зависимым признаком методом попарного удаления "NA"
b <- as.data.frame(na.omit(c))
groups <- data.frame(table(b$factor)) #Определяем "n" и переводим таблицу в объект "data.frame"
names(groups) <- c("id","n") #Переименовываем созданную таблицу - "id","n"
groups <- subset(groups, n >= ingroup, select = c(id,n)) #Выбираем только группы, в которых не менее ingroup вариант
ifgroups <- c$factor %in% groups$id #Создаем новый вектор "ifgroups", где TRUE >= ingroup и FALSE < ingroup
lgroups <- as.data.frame(subset(c, ifgroups == "TRUE"))
return(lgroups) #Возвращаем таблицу с нужной величиной градации
})}
# Источник: Hozo, S.P. Estimating the mean and variance from the median, range, and the size of a sample
# / B. Djulbegovic & I. Hozo // BMC Medical Research Methodology. – 2005. – Vol. 5, – nr 1. – P. 13.
mean.sx2 <- function(x)
{
a <- min(x)
b <- max(x)
m <- median(x)
n <- sum(!is.na(x))
mn <- (a+2*m+b)/4+(a-2*m+b)/(4*n)
s <- sqrt((a^2+m^2+b^2+(n-3)*((a+m)^2+(m+b)^2)/8-n*mn^2)/(n-1))
sx <- s/sqrt(sum(!is.na(x)))
c(mn,sx)
}
#Автоматическое округление в зависимости от исходных данных
round_auto <- function(x){
ifelse(x < 1, x <- round(x, 3),
ifelse(nchar(trunc(x)) >= 4, x <- round(x),
ifelse(nchar(trunc(x)) < 2, x <- round(x, 2), x <- round(x, 1))
)
)
return(x) #На выходе - округлённое значение
}
se <- function(x)
{
n <- sum(!is.na(x))
se_output <- sd(na.omit(x))/sqrt(n)
return(se_output)
}
library(nortest)
#Here Mean_ng, SE_ng mean nongaussian distribution
descrstats2 <- function(x, grx, nmin = 3, Mean = TRUE, SE = TRUE, Mean_ng = TRUE, SE_ng = TRUE, Me = TRUE, Min = TRUE, Max = TRUE,
Range = TRUE, As = FALSE, Ex = FALSE, Q1 = TRUE, Q3 = TRUE, IQR = TRUE, SD = TRUE,
Cv = TRUE, AD = FALSE, SF = FALSE, srclevs = FALSE)
{
descrtable<-function(x = x, grx = grx, nmin = nmin, Mean_ = Mean, SE_ = SE, Mean_ng_ = Mean_ng, SE_ng_ = SE_ng, Me_ = Me,
Min_ = Min, Max_ = Max, Range_ = Range, As_ = As, Ex_ = Ex, Q1_ = Q1, Q3_ = Q3, IQR_ = IQR, SD_ = SD,
Cv_ = Cv, AD_ = AD, SF_ = SF, srclevs = srclevs)
{
grx <- as.data.frame(interaction(grx, sep = ":"))
xname <- names(x)
nminame <- names(grx)
x <- as.data.frame(x)
x <- cbind(x,grx)
#Удаляем строки, где количество вариант меньше определенного grf минимума
#x - таблица: зав.пр, групп.пр; кол-во вариант в группе - nmin
x <- n_group_restricted(x, nmin)
names(x) <- c("x","grx")
x <- na.omit(x)
res <- round_auto(as.data.frame(tapply(x$x, x$grx, FUN=function(y) sum(!is.na(y)))))
names(res) <- "n"
#res$naa<-apply(res$n,1,is.na)
if(Mean_) {res$Mean <- tapply(x$x, x$grx, FUN = function(y) round_auto(mean(y, na.rm=TRUE)))}
if(SE_) {res$SE <- tapply(x$x, x$grx, FUN=function(y) round_auto(se(y)))}
if(Mean_ng_) {res$Mean_ng <- tapply(x$x, x$grx, FUN = function(y) round_auto(mean.sx2(y)[1]))}
if(SE_ng_) {res$SE_ng <- tapply(x$x, x$grx, FUN = function(y) round_auto(mean.sx2(y)[2]))}
if(Me_) {res$Me<-tapply(x$x,x$grx,FUN=function(y) round_auto(median(y)))}
if(Min_) {res$Min<-tapply(x$x,x$grx,FUN=function(y) round_auto(min(y)))}
if(Max_) {res$Max<-tapply(x$x,x$grx,FUN=function(y) round_auto(max(y)))}
if(Range_) {res$Range<-tapply(x$x,x$grx,FUN=function(y) round_auto(max(y)-min(y)))}
if(As_) {res$As<-tapply(x$x,x$grx,FUN=function(y) round_auto(skewness(y)))}
if(Ex_) {res$Ex<-tapply(x$x,x$grx,FUN=function(y) round_auto(kurtosis(y)))}
if(Q1_) {res$Q1<-tapply(x$x,x$grx,FUN=function(y) round_auto(quantile(y, probs=.25,na.rm = TRUE, type = 8)))}
if(Q3_) {res$Q3<-tapply(x$x,x$grx,FUN=function(y) round_auto(quantile(y, probs=.75,na.rm = TRUE, type = 8)))}
if(IQR_) {res$IQR<-tapply(x$x,x$grx,FUN=function(y)
round_auto(quantile(y, probs=.75,na.rm = TRUE, type = 8)-quantile(y, probs=.25,na.rm = TRUE, type = 8)))}
if(SD_) {res$SD<-tapply(x$x,x$grx,FUN=function(y) round_auto(sd(y)))}
if(Cv_) {res$Cv<-tapply(x$x,x$grx,FUN=function(y) round((sd(y)*100/mean(y)), 1))}
if(AD_) {res$AD<-tapply(x$x,x$grx,FUN=function(y)
ifelse(is.na(sum(!is.na(y))),NA, ifelse(sum(!is.na(y))>7, round(ad.test(y)$statistic,3),NA)))
res$AD.p<-tapply(x$x,x$grx,FUN=function(y)
ifelse(is.na(sum(!is.na(y))),NA, ifelse(sum(!is.na(y))>7, round(ad.test(y)$p,3),NA)))}
if(SF_) {res$SF<-tapply(x$x,x$grx,FUN=function(y)
ifelse(is.na(sum(!is.na(y))),NA, ifelse(sum(!is.na(y))>7, round(sf.test(y)$statistic,3),NA)))
res$SF.p<-tapply(x$x,x$grx,FUN=function(y)
ifelse(is.na(sum(!is.na(y))),NA, ifelse(sum(!is.na(y))>7, round(sf.test(y)$p,3),NA)))}
res<-cbind(Group=row.names(res),res)
rn<-paste(xname,1:nrow(res),sep="_")
# row.names(res)<-rn
Index<-rep(xname,nrow(res))
res<-data.frame(Index,res)#Добавляем стобец с повторяющимися названиями строк
res<-subset(res,is.na(res$n)==FALSE)
names(res)<-c("Index",names(res[2:ncol(res)]))
row.names(res)<-NULL
return(res)
}
if(ncol(x)!=1){
res<-descrtable(x[1],grx,nmin)
for(i in 2:ncol(x))
{
res2<-descrtable(x[i],grx,nmin)
res<-rbind(res,res2)
}
return(res)
}
}Пример использования
iris$factor <- gl(3, 1, 150)
descrstats2(iris[1:4], list(iris[,5], iris[,6]), Mean_ng = FALSE, SE_ng = FALSE, Q1 = FALSE, Q3 = FALSE, Me = FALSE)## Index Group n Mean SE Min Max Range IQR SD Cv
## 1 Sepal.Length setosa:1 17 5.050 0.095 4.4 5.7 1.3 0.567 0.391 7.7
## 2 Sepal.Length versicolor:1 17 5.770 0.131 4.9 6.6 1.7 0.833 0.538 9.3
## 3 Sepal.Length virginica:1 16 6.760 0.147 5.8 7.7 1.9 0.817 0.588 8.7
## 4 Sepal.Length setosa:2 17 5.010 0.067 4.3 5.4 1.1 0.233 0.276 5.5
## 5 Sepal.Length versicolor:2 16 6.020 0.105 5.6 6.9 1.3 0.458 0.418 6.9
## 6 Sepal.Length virginica:2 17 6.450 0.164 4.9 7.7 2.8 0.667 0.676 10.5
## 7 Sepal.Length setosa:3 16 4.950 0.099 4.4 5.8 1.4 0.517 0.395 8.0
## 8 Sepal.Length versicolor:3 17 6.020 0.137 5.1 7.0 1.9 1.000 0.564 9.4
## 9 Sepal.Length virginica:3 17 6.570 0.155 5.7 7.9 2.2 0.867 0.639 9.7
## 10 Sepal.Width setosa:1 17 3.470 0.097 3.0 4.4 1.4 0.600 0.400 11.5
## 11 Sepal.Width versicolor:1 17 2.680 0.080 2.0 3.2 1.2 0.533 0.329 12.3
## 12 Sepal.Width virginica:1 16 2.980 0.074 2.5 3.8 1.3 0.258 0.297 10.0
## 13 Sepal.Width setosa:2 17 3.450 0.070 3.0 3.9 0.9 0.400 0.290 8.4
## 14 Sepal.Width versicolor:2 16 2.910 0.052 2.6 3.4 0.8 0.258 0.208 7.2
## 15 Sepal.Width virginica:2 17 3.020 0.075 2.5 3.6 1.1 0.500 0.309 10.2
## 16 Sepal.Width setosa:3 16 3.360 0.112 2.3 4.1 1.8 0.475 0.450 13.4
## 17 Sepal.Width versicolor:3 17 2.740 0.085 2.2 3.3 1.1 0.567 0.352 12.9
## 18 Sepal.Width virginica:3 17 2.920 0.089 2.2 3.8 1.6 0.467 0.366 12.5
## 19 Petal.Length setosa:1 17 1.490 0.035 1.3 1.9 0.6 0.100 0.145 9.8
## 20 Petal.Length versicolor:1 17 4.200 0.119 3.3 4.9 1.6 0.667 0.490 11.7
## 21 Petal.Length virginica:1 16 5.570 0.147 4.8 6.7 1.9 0.758 0.586 10.5
## 22 Petal.Length setosa:2 17 1.410 0.040 1.0 1.6 0.6 0.133 0.165 11.7
## 23 Petal.Length versicolor:2 16 4.330 0.105 3.5 4.9 1.4 0.558 0.421 9.7
## 24 Petal.Length virginica:2 17 5.490 0.138 4.5 6.9 2.4 0.700 0.569 10.4
## 25 Petal.Length setosa:3 16 1.490 0.052 1.2 1.9 0.7 0.358 0.206 13.9
## 26 Petal.Length versicolor:3 17 4.250 0.124 3.0 5.1 2.1 0.700 0.511 12.0
## 27 Petal.Length virginica:3 17 5.600 0.128 5.0 6.7 1.7 0.833 0.529 9.4
## 28 Petal.Width setosa:1 17 0.229 0.021 0.1 0.4 0.3 0.100 0.085 37.0
## 29 Petal.Width versicolor:1 17 1.290 0.050 1.0 1.5 0.5 0.433 0.205 15.8
## 30 Petal.Width virginica:1 16 2.070 0.066 1.6 2.5 0.9 0.500 0.263 12.7
## 31 Petal.Width setosa:2 17 0.247 0.030 0.1 0.6 0.5 0.100 0.123 49.8
## 32 Petal.Width versicolor:2 16 1.340 0.051 1.0 1.8 0.8 0.217 0.203 15.2
## 33 Petal.Width virginica:2 17 2.090 0.071 1.5 2.5 1.0 0.433 0.291 13.9
## 34 Petal.Width setosa:3 16 0.262 0.027 0.1 0.5 0.4 0.158 0.109 41.4
## 35 Petal.Width versicolor:3 17 1.350 0.047 1.0 1.7 0.7 0.300 0.194 14.4
## 36 Petal.Width virginica:3 17 1.920 0.062 1.4 2.4 1.0 0.233 0.254 13.2
No comments:
Post a Comment