[81de4e]: / weighted.stats.R

Download this file

99 lines (79 with data), 2.9 kB

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
weighted.stats =
function (x, w, conditions)
{
n = ncol(x) #number of samples
p = nrow(x) #number of genes
nclass = length(unique(conditions)) #number of class
if(is.factor(conditions)) {cNames = sort(levels(conditions))}
if (is.numeric(conditions)) {cNames = as.character(sort(unique(conditions)))}
WM = WS = wSum = se.scale = matrix(0, p, nclass)
rownames(WS) = rownames(WM) = rownames(wSum) = rownames(se.scale) = rownames(x)
colnames(WS) = colnames(WM) = colnames(wSum) = colnames(se.scale) = cNames
c.ind = as.numeric(conditions)
w.mean00 =
function (x, w)
{
wm = NULL
for (i in 1:p)
{
wm0 = sum(w[i,]*x[i,]) / sum(w[i,])
wm = c(wm, wm0)
}
return(wm)
}
w.mean =
function (x, w, conditions)
{
for (j in 1:nclass)
{
WM[,j] = w.mean00(x[,c.ind == j], w[,c.ind == j])
}
return(WM)
}
w.sd =
function (x, w, conditions)
{
w.sd00 =
function (x, w)
{
ws = NULL
w.sd0 =
function (x, w)
{
sumw = sum(w)
sumw.sq = sum(w)^2
w.sq = sum(w^2)
denom = sum(w * ((x - mean(x))^2))
sqrt((sumw * denom) / (sumw.sq - w.sq))
}
for (i in 1:p)
{
ws0 = w.sd0(x[i,], w[i,])
ws = c(ws, ws0)
}
return(ws)
}
for (j in 1:nclass)
{
WS[,j] = w.sd00(x[,c.ind == j], w[,c.ind == j])
}
return(WS)
}
weightedMean = w.mean00(x, w) #Overall weighted mean
weightedMean.C = w.mean(x, w, conditions) #Weighted means for each group
weightedSD.C = w.sd(x, w, conditions) #Weighted standard deviations for each group
#weightedSD.pooled = weightedSD.C
for (i in 1:nclass)
{
tmp = w[,which(c.ind == i)]
rSum.tmp = rowSums(tmp)
wSum[,i] = rSum.tmp
}
weightedSD.pooled = sqrt(rowSums((wSum-1) * (weightedSD.C^2)) / (rowSums(wSum) - nclass))
se.scale = sqrt(1 / wSum + 1 / rowSums(wSum))
s0 = median(weightedSD.pooled)
delta = (weightedMean.C - weightedMean)/(se.scale*(weightedSD.pooled + s0))
weightedSD.pooled = sqrt(rowSums(as.data.frame(weightedSD.pooled)) / (n - nclass))
stats = list(n = n, p = p, nclass = nclass, se.scale = se.scale, weightedMean = weightedMean, weightedMean.C = weightedMean.C, weightedSD.C = weightedSD.C, weightedSD.pooled = weightedSD.pooled, delta = delta)
return(stats)
}