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