関数を作ってみました。
HS.tapply <- function(x, index, type) { # create index indicate matrix index <- as.numeric(factor(index)) mindex <- NULL for(i in 1:length(unique(index))) { mindex <- cbind(mindex, index-(i-1)) } if (type == 5 || type == "max" || type == 6 || type == "min") { mindex[mindex!=1] <- NA } else { mindex[mindex!=1] <- 0 } r <- NULL # Mean if (type == 1 || type == "mean") { r <- t(mindex) %*% x / as.numeric(table(index)) } # Summation else if (type == 2 || type == "sum") { r <- t(mindex) %*% x } # Sum of squared deviation else if (type == 3 || type == "ss") { r <- t(mindex) %*% x^2 - (t(mindex) %*% x)^2 / as.numeric(table(index)) } # Variance (sample variance) else if (type == 4 || type == "var") { r <- (t(mindex) %*% x^2 - (t(mindex) %*% x)^2 / as.numeric(table(index))) / (as.numeric(table(index))-1) } # Maximum else if (type == 5 || type == "max") { for(i in 1:length(unique(index))) { r <- rbind(r, x[cbind( max.col(ifelse(is.na(t(mindex[,i]*x)), -Inf, t(mindex[,i]*x))), 1:(ncol(x)) )]) } } # Minimum else if (type == 6 || type == "min") { for(i in 1:length(unique(index))) { r <- rbind(r, x[cbind( max.col(ifelse(is.na(t(mindex[,i]*x)), -Inf, -t(mindex[,i]*x))), 1:(ncol(x)) )]) } } return(r) }
データが行列状になっていても計算してくれます(グループ変数が等しく、異なる測定値がたくさんあるようなデータ)。
データが多いとメモリが足りなくなるので、適当に分割して計算します。
n <- c(5, 4, 7) g <- rep(1:length(n), n) gs <- c("a", "a", "a", "a", "a", "b", "b", "b", "b", "c", "c", "c", "c", "c", "c", "c") rep <- 20000 x <- matrix(rnorm(sum(n)*rep), ncol=rep, nrow=sum(n)) Rprof() varA <- apply(x, 2, function(x, g) tapply(x, g, var), g) Rprof(NULL) print(a<-summaryRprof()$sampling.time) Rprof() varF <- HS.tapply(x, g, "var") #varF <- HS.tapply(x, gs, "var") # グループ変数は文字列等でもOK Rprof(NULL) print(f<-summaryRprof()$sampling.time) max(varA-varF) varA[,1:3] varF[,1:3] f/a
あまりの早さにびっくり。
> print(a<-summaryRprof()$sampling.time) [1] 12.02 > print(f<-summaryRprof()$sampling.time) [1] 0.02 > max(varA-varF) [1] 1.776357e-15 > varA[,1:3] [,1] [,2] [,3] 1 0.70445164 2.2072017 0.6809683 2 0.09444263 0.8143882 0.4464633 3 0.28558027 0.9014923 0.9471029 > varF[,1:3] [,1] [,2] [,3] [1,] 0.70445164 2.2072017 0.6809683 [2,] 0.09444263 0.8143882 0.4464633 [3,] 0.28558027 0.9014923 0.9471029 > f/a [1] 0.001663894
条件演算子
R の条件演算子には &, |, &&, || などがありますが、前者二つはベクトル化されている演算子であり、後者二つはベクトル化されていないという違いがあります。
したがって、データフレームなどに subset 関数をつかって条件を満たす症例だけを抽出したい場合などは両者を混同しないように注意が必要です。
df<-data.frame(x=runif(10),y=rnorm(10)) subset(df, y<0 & x>0.5) # (a) & を使う subset(df, y<0 && x>0.5) # (b) && を使う
(a) は正しく二つの条件を同時に満たす症例を抽出できますが、(b) はそうなりません。
