Browse code

add R level interfaces for fitting a row-col effects model (ala RMA) to a user specified matrix

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/preprocessCore@27176 bc3139a8-67e5-0310-9ffc-ced21a209358

Ben Bolstad authored on 15/09/2007 23:27:23
Showing6 changed files

... ...
@@ -1,12 +1,12 @@
1 1
 Package: preprocessCore
2
-Version: 0.99.16
2
+Version: 0.99.17
3 3
 Title: A collection of pre-processing functions
4 4
 Author: Benjamin Milo Bolstad <bmb@bmbolstad.com>
5 5
 Maintainer: Benjamin Milo Bolstad <bmb@bmbolstad.com>
6 6
 Depends: methods
7 7
 Description: A library of core preprocessing routines 
8 8
 License: LGPL version 2 or newer
9
-Collate:  normalize.quantiles.R quantile_extensions.R init.R
9
+Collate:  normalize.quantiles.R quantile_extensions.R rcModel.R init.R
10 10
 LazyLoad: yes
11 11
 biocViews: Infrastructure
12 12
 
13 13
new file mode 100644
... ...
@@ -0,0 +1,47 @@
1
+
2
+
3
+
4
+rcModelPLM <- function(y){
5
+  if (!is.matrix(y))
6
+    stop("argument should be matrix")
7
+  PsiCode <- 0
8
+  PsiK <- 1.345
9
+  .Call("R_rlm_rma_default_model",y,PsiCode,PsiK,PACKAGE="preprocessCore")
10
+}
11
+
12
+
13
+
14
+rcModelWPLM <- function(y, w){
15
+  if (!is.matrix(y))
16
+    stop("argument should be matrix")
17
+  if (is.vector(w)){
18
+    if (length(w) != prod(dim(y))){
19
+      stop("weights are not correct length")
20
+    }
21
+  } else if (is.matrix(w)){
22
+    if (!all(dim(w) == dim(y))){
23
+      stop("weights should be same dimension as input matrix")
24
+    }
25
+
26
+  }
27
+  if (any(w < 0)){
28
+    stop("weights should be no negative")
29
+  }
30
+
31
+  
32
+    
33
+  PsiCode <- 0
34
+  PsiK <- 1.345
35
+  .Call("R_wrlm_rma_default_model",y,PsiCode,PsiK,as.double(w),PACKAGE="preprocessCore")
36
+
37
+}
38
+
39
+
40
+
41
+rcModelMedianPolish <- function(y){
42
+  if (!is.matrix(y))
43
+    stop("argument should be matrix")
44
+  PsiCode <- 0
45
+  PsiK <- 1.345
46
+  .Call("R_medianpolish_rma_default_model",y,PACKAGE="preprocessCore")
47
+}
0 48
new file mode 100644
... ...
@@ -0,0 +1,62 @@
1
+\name{rcModels}
2
+\alias{rcModelPLM}
3
+\alias{rcModelWPLM}
4
+\alias{rcModelMedianPolish}
5
+\title{Fit row-column model to a matrix}
6
+\description{These functions fit row-column effect models to matrices
7
+}
8
+\usage{
9
+rcModelPLM(y)
10
+rcModelWPLM(y, w)
11
+rcModelMedianPolish(y)
12
+}
13
+\arguments{
14
+  \item{y}{A matrix}
15
+  \item{w}{A matrix or vector of weights. These should be non-negative.}
16
+}
17
+\value{
18
+  A list with following items:
19
+  \item{Estimates}{The parameter estimates. Stored in column effect then
20
+    row effect order}
21
+  \item{Weights}{The final weights used}
22
+  \item{Residuals}{The residuals}
23
+  \item{StdErrors}{Standard error estimates. Stored in column effect
24
+    then row effect order}
25
+}
26
+\details{
27
+  These functions fit row-column models to the specified input
28
+  matrix. Specifically the model \deqn{y_{ij} = r_i + c_j +
29
+    \epsilon_{ij}}{y_ij = r_i + c_j + e_ij}
30
+       with \eqn{r_i} and \eqn{c_j} as row and column effects
31
+       respectively. Note that this functions treat the row effect as
32
+       the parameter to be constrained using sum to zero (for
33
+       \code{rcModelPLM} and \code{rcModelWPLM}) or median of zero (for
34
+       \code{rcModelMedianPolish}).
35
+
36
+       The \code{rcModelPLM} and \code{rcModelWPLM} functions use a
37
+       robust linear model procedure for fitting the model.
38
+
39
+       The function \code{rcModelMedianPolish} uses the median polish algorithm.
40
+}
41
+\seealso{
42
+}
43
+\examples{
44
+col.effects <- c(10,11,10.5,12,9.5)
45
+row.effects <- c(seq(-0.5,-0.1,by=0.1),seq(0.1,0.5,by=0.1))
46
+
47
+
48
+y <- outer(row.effects, col.effects,"+")
49
+w <- runif(50)
50
+
51
+rcModelPLM(y)
52
+rcModelWPLM(y, w)
53
+rcModelMedianPolish(y)
54
+
55
+y <- y + rnorm(50)
56
+
57
+rcModelPLM(y)
58
+rcModelWPLM(y, w)
59
+rcModelMedianPolish(y)
60
+}
61
+\author{B. M. Bolstad \email{bmb@bmbolstad.com}}
62
+\keyword{models}
0 63
\ No newline at end of file
... ...
@@ -4,7 +4,7 @@
4 4
  **
5 5
  ** Aim: Code which provides .Call() interfaces to the rlm code.
6 6
  **
7
- ** Copyright (C) 2006 Ben Bolstad
7
+ ** Copyright (C) 2006-2007 Ben Bolstad
8 8
  **
9 9
  ** created by: B. M. Bolstad <bmb@bmbolstad.com>
10 10
  ** 
... ...
@@ -13,7 +13,9 @@
13 13
  ** History
14 14
  ** Aug 16, 2006 Initial version
15 15
  ** Nov 1, 2006 - add SE to output of function
16
- **
16
+ ** Sep 13, 2007 - Make the value of the constrained parameters something more sensible
17
+ ** Sep 14, 2007 - Add medianpolish code interface (yes it is not really an rlm method, 
18
+ **                but it is analogous enough in the format presented here)
17 19
  **
18 20
  **
19 21
  *********************************************************************/
... ...
@@ -27,6 +29,9 @@
27 29
 
28 30
 #include "psi_fns.h"
29 31
 
32
+#include "medianpolish.h"
33
+
34
+
30 35
 #include "R_rlm_interfaces.h"
31 36
 
32 37
 
... ...
@@ -74,6 +79,8 @@ SEXP R_rlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK){
74 79
 
75 80
   int rows;
76 81
   int cols;
82
+
83
+  int i;
77 84
   
78 85
   PROTECT(dim1 = getAttrib(Y,R_DimSymbol));
79 86
   rows = INTEGER(dim1)[0];
... ...
@@ -111,8 +118,13 @@ SEXP R_rlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK){
111 118
 
112 119
   
113 120
 
121
+  beta[rows+cols -1] = 0.0;
122
+  se[rows+cols -1] = 0.0;
114 123
 
115 124
 
125
+  for (i = cols; i < rows + cols -1; i++)
126
+    beta[rows+cols -1]-=beta[i];
127
+
116 128
 
117 129
 
118 130
 
... ...
@@ -165,6 +177,8 @@ SEXP R_wrlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK, SEXP Weights){
165 177
 
166 178
   int rows;
167 179
   int cols;
180
+
181
+  int i;
168 182
   
169 183
   PROTECT(dim1 = getAttrib(Y,R_DimSymbol));
170 184
   rows = INTEGER(dim1)[0];
... ...
@@ -196,16 +210,120 @@ SEXP R_wrlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK, SEXP Weights){
196 210
 
197 211
   
198 212
   rlm_wfit_anova(Ymat, rows, cols, w, beta, residuals, weights, PsiFunc(asInteger(PsiCode)),asReal(PsiK), 20, 0);
199
-  
200 213
   rlm_compute_se_anova(Ymat, rows, cols, beta, residuals, weights,se, (double *)NULL, &residSE, 4, PsiFunc(asInteger(PsiCode)),asReal(PsiK));
201 214
 
202 215
 
216
+  beta[rows+cols -1] = 0.0;
217
+  se[rows+cols -1] = 0.0;
218
+
219
+  for (i = cols; i < rows + cols -1; i++)
220
+    beta[rows+cols -1]-=beta[i];
221
+
222
+
223
+
224
+
225
+
226
+
227
+
228
+  PROTECT(R_return_value_names= allocVector(STRSXP,4));
229
+  SET_VECTOR_ELT(R_return_value_names,0,mkChar("Estimates"));
230
+  SET_VECTOR_ELT(R_return_value_names,1,mkChar("Weights"));
231
+  SET_VECTOR_ELT(R_return_value_names,2,mkChar("Residuals"));
232
+  SET_VECTOR_ELT(R_return_value_names,3,mkChar("StdErrors"));
233
+  setAttrib(R_return_value, R_NamesSymbol,R_return_value_names);
234
+  UNPROTECT(2);
235
+  return R_return_value;
236
+
237
+}
238
+
239
+
240
+
241
+
242
+
243
+
244
+
245
+
246
+
247
+
248
+
249
+
250
+
251
+
252
+
253
+
254
+
255
+
256
+
257
+SEXP R_medianpolish_rma_default_model(SEXP Y){
258
+
259
+
260
+
261
+  SEXP R_return_value;
262
+  SEXP R_weights;
263
+  SEXP R_residuals;
264
+  SEXP R_beta;
265
+  SEXP R_SE;
203 266
   
267
+  SEXP R_return_value_names;
204 268
 
269
+  SEXP dim1;
205 270
 
271
+  double *beta;
272
+  double *residuals;
273
+  double *weights;
274
+  double *se;
206 275
 
276
+  double intercept;
207 277
 
278
+  double *Ymat;
208 279
 
280
+  int rows;
281
+  int cols;
282
+
283
+  int i;
284
+  
285
+  PROTECT(dim1 = getAttrib(Y,R_DimSymbol));
286
+  rows = INTEGER(dim1)[0];
287
+  cols = INTEGER(dim1)[1];
288
+  UNPROTECT(1);
289
+
290
+  PROTECT(R_return_value = allocVector(VECSXP,4));
291
+  PROTECT(R_beta = allocVector(REALSXP, rows + cols));
292
+  /* PROTECT(R_weights = allocMatrix(REALSXP,rows,cols));*/
293
+  PROTECT(R_residuals = allocMatrix(REALSXP,rows,cols));
294
+  /*  PROTECT(R_SE = allocVector(REALSXP,rows+cols)); */
295
+
296
+  R_weights = R_NilValue;
297
+  R_SE = R_NilValue;
298
+
299
+
300
+  SET_VECTOR_ELT(R_return_value,0,R_beta);
301
+  SET_VECTOR_ELT(R_return_value,1,R_weights);
302
+  SET_VECTOR_ELT(R_return_value,2,R_residuals);
303
+  SET_VECTOR_ELT(R_return_value,3,R_SE);
304
+
305
+  UNPROTECT(2);
306
+
307
+  beta = NUMERIC_POINTER(R_beta);
308
+  residuals = NUMERIC_POINTER(R_residuals);
309
+  /*  weights = NUMERIC_POINTER(R_weights);
310
+      se = NUMERIC_POINTER(R_SE);
311
+  */
312
+
313
+  Ymat = NUMERIC_POINTER(Y);
314
+
315
+  for (i=0; i < rows*cols; i++){
316
+    residuals[i] = Ymat[i];
317
+  }
318
+
319
+  memset(beta, 0, (rows+cols)*sizeof(double));
320
+
321
+  
322
+  
323
+  median_polish_fit_no_copy(residuals, rows, cols, &beta[cols], &beta[0], &intercept);
324
+
325
+  for (i=0; i < cols; i++)
326
+    beta[i]+=intercept;
209 327
 
210 328
 
211 329
 
... ...
@@ -219,4 +337,12 @@ SEXP R_wrlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK, SEXP Weights){
219 337
   UNPROTECT(2);
220 338
   return R_return_value;
221 339
 
340
+
341
+
342
+
343
+
344
+
345
+
346
+
347
+
222 348
 }
... ...
@@ -4,5 +4,6 @@
4 4
 
5 5
 SEXP R_rlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK);
6 6
 SEXP R_wrlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK, SEXP Weights);
7
+SEXP R_medianpolish_rma_default_model(SEXP Y);
7 8
 
8 9
 #endif
... ...
@@ -49,6 +49,7 @@ static const R_CallMethodDef callMethods[]  = {
49 49
   {"R_qnorm_within_blocks",(DL_FUNC)&R_qnorm_within_blocks,3},
50 50
   {"R_rlm_rma_default_model",(DL_FUNC)&R_rlm_rma_default_model,3},
51 51
   {"R_wrlm_rma_default_model", (DL_FUNC)&R_wrlm_rma_default_model,4},
52
+  {"R_medianpolish_rma_default_model", (DL_FUNC)&R_medianpolish_rma_default_model,1},
52 53
   {NULL, NULL, 0}
53 54
   };
54 55