Ben Bolstad
| Snapshot Date: 2020-10-16 14:40:19 -0400 (Fri, 16 Oct 2020) |
| URL: https://git.bioconductor.org/packages/affyPLM |
| Branch: RELEASE_3_11 |
| Last Commit: 11b83c7 |
| Last Changed Date: 2020-04-27 14:08:49 -0400 (Mon, 27 Apr 2020) |
| Back to Multiple platform build/check report for BioC 3.11 |
|
This page was generated on 2020-10-17 11:54:21 -0400 (Sat, 17 Oct 2020).
| TO THE DEVELOPERS/MAINTAINERS OF THE affyPLM PACKAGE: Please make sure to use the following settings in order to reproduce any error or warning you see on this page. |
| Package 34/1905 | Hostname | OS / Arch | INSTALL | BUILD | CHECK | BUILD BIN | ||||||
| affyPLM 1.64.0 Ben Bolstad
| malbec2 | Linux (Ubuntu 18.04.4 LTS) / x86_64 | OK | OK | [ OK ] | | ||||||
| tokay2 | Windows Server 2012 R2 Standard / x64 | OK | OK | OK | OK | | ||||||
| machv2 | macOS 10.14.6 Mojave / x86_64 | OK | OK | OK | OK | |
| Package: affyPLM |
| Version: 1.64.0 |
| Command: /home/biocbuild/bbs-3.11-bioc/R/bin/R CMD check --install=check:affyPLM.install-out.txt --library=/home/biocbuild/bbs-3.11-bioc/R/library --no-vignettes --timings affyPLM_1.64.0.tar.gz |
| StartedAt: 2020-10-16 22:42:43 -0400 (Fri, 16 Oct 2020) |
| EndedAt: 2020-10-16 22:46:20 -0400 (Fri, 16 Oct 2020) |
| EllapsedTime: 217.4 seconds |
| RetCode: 0 |
| Status: OK |
| CheckDir: affyPLM.Rcheck |
| Warnings: 0 |
##############################################################################
##############################################################################
###
### Running command:
###
### /home/biocbuild/bbs-3.11-bioc/R/bin/R CMD check --install=check:affyPLM.install-out.txt --library=/home/biocbuild/bbs-3.11-bioc/R/library --no-vignettes --timings affyPLM_1.64.0.tar.gz
###
##############################################################################
##############################################################################
* using log directory ‘/home/biocbuild/bbs-3.11-bioc/meat/affyPLM.Rcheck’
* using R version 4.0.3 (2020-10-10)
* using platform: x86_64-pc-linux-gnu (64-bit)
* using session charset: UTF-8
* using option ‘--no-vignettes’
* checking for file ‘affyPLM/DESCRIPTION’ ... OK
* this is package ‘affyPLM’ version ‘1.64.0’
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package ‘affyPLM’ can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking ‘build’ directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking use of PKG_*FLAGS in Makefiles ... OK
* checking compiled code ... NOTE
Note: information on .o files is not available
* checking files in ‘vignettes’ ... OK
* checking examples ... OK
Examples with CPU (user + system) or elapsed time > 5s
user system elapsed
threestep 18.620 0.007 18.631
fitPLM 12.353 0.155 12.508
PLMset2exprSet 6.726 0.119 6.847
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ...
Running ‘C_code_tests.R’
Running ‘PLM_tests.R’
Running ‘preprocess_tests.R’
Running ‘threestepPLM_tests.R’
OK
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in ‘inst/doc’ ... OK
* checking running R code from vignettes ... SKIPPED
* checking re-building of vignette outputs ... SKIPPED
* checking PDF version of manual ... OK
* DONE
Status: 1 NOTE
See
‘/home/biocbuild/bbs-3.11-bioc/meat/affyPLM.Rcheck/00check.log’
for details.
affyPLM.Rcheck/00install.out
##############################################################################
##############################################################################
###
### Running command:
###
### /home/biocbuild/bbs-3.11-bioc/R/bin/R CMD INSTALL affyPLM
###
##############################################################################
##############################################################################
* installing to library ‘/home/biocbuild/bbs-3.11-bioc/R/library’
* installing *source* package ‘affyPLM’ ...
** using staged installation
** libs
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c LESN.c -o LESN.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_avg_log.c -o PLM_avg_log.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_biweight.c -o PLM_biweight.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_log_avg.c -o PLM_log_avg.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_medianPM.c -o PLM_medianPM.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_median_logPM.c -o PLM_median_logPM.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_medianpolish.c -o PLM_medianpolish.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c PLM_modelmatrix.c -o PLM_modelmatrix.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c SCAB.c -o SCAB.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c chipbackground.c -o chipbackground.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c common_types.c -o common_types.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c do_PLMrlm.c -o do_PLMrlm.o
do_PLMrlm.c: In function ‘do_PLM_rlm’:
do_PLMrlm.c:620:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c do_PLMrma.c -o do_PLMrma.o
do_PLMrma.c: In function ‘do_PLMrma’:
do_PLMrma.c:209:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c do_PLMthreestep.c -o do_PLMthreestep.o
do_PLMthreestep.c: In function ‘do_PLMthreestep’:
do_PLMthreestep.c:118:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c idealmismatch.c -o idealmismatch.o
idealmismatch.c: In function ‘IdealMM_correct_single’:
idealmismatch.c:71:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
idealmismatch.c: In function ‘SpecificBiweightCorrect_single’:
idealmismatch.c:183:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c lm_threestep.c -o lm_threestep.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c matrix_functions.c -o matrix_functions.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c nthLargestPM.c -o nthLargestPM.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c preprocess.c -o preprocess.o
preprocess.c: In function ‘pp_background’:
preprocess.c:158:7: warning: variable ‘which_lesn’ set but not used [-Wunused-but-set-variable]
int which_lesn;
^~~~~~~~~~
preprocess.c: In function ‘pp_bothstages’:
preprocess.c:677:12: warning: variable ‘cols’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
preprocess.c:677:7: warning: variable ‘rows’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c psi_fns.c -o psi_fns.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c qnorm_probeset.c -o qnorm_probeset.o
qnorm_probeset.c: In function ‘qnorm_probeset_c’:
qnorm_probeset.c:110:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c read_rmaexpress.c -o read_rmaexpress.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c rlm_PLM.c -o rlm_PLM.o
rlm_PLM.c: In function ‘R_rlm_PLMset_c’:
rlm_PLM.c:1481:12: warning: variable ‘cols’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
rlm_PLM.c:1481:7: warning: variable ‘rows’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c rlm_threestep.c -o rlm_threestep.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c rmaPLM_pseudo.c -o rmaPLM_pseudo.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c rma_PLM.c -o rma_PLM.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c rma_common.c -o rma_common.o
rma_common.c: In function ‘median’:
rma_common.c:60:7: warning: unused variable ‘i’ [-Wunused-variable]
int i;
^
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c scaling.c -o scaling.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c threestep.c -o threestep.o
threestep.c: In function ‘threestep_summary’:
threestep.c:82:15: warning: variable ‘MM’ set but not used [-Wunused-but-set-variable]
double *PM,*MM;
^~
threestep.c: In function ‘R_threestep_c’:
threestep.c:193:12: warning: variable ‘cols’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
threestep.c:193:7: warning: variable ‘rows’ set but not used [-Wunused-but-set-variable]
int rows,cols;
^~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c threestep_PLM.c -o threestep_PLM.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c threestep_common.c -o threestep_common.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c threestep_summary.c -o threestep_summary.o
threestep_summary.c: In function ‘do_3summary’:
threestep_summary.c:73:7: warning: variable ‘first_ind’ set but not used [-Wunused-but-set-variable]
int first_ind;
^~~~~~~~~
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c threestep_summary_methods.c -o threestep_summary_methods.o
gcc -I"/home/biocbuild/bbs-3.11-bioc/R/include" -DNDEBUG -I'/home/biocbuild/bbs-3.11-bioc/R/library/preprocessCore/include' -I/usr/local/include -fpic -g -O2 -Wall -c transfns.c -o transfns.o
gcc -shared -L/home/biocbuild/bbs-3.11-bioc/R/lib -L/usr/local/lib -o affyPLM.so LESN.o PLM_avg_log.o PLM_biweight.o PLM_log_avg.o PLM_medianPM.o PLM_median_logPM.o PLM_medianpolish.o PLM_modelmatrix.o SCAB.o chipbackground.o common_types.o do_PLMrlm.o do_PLMrma.o do_PLMthreestep.o idealmismatch.o lm_threestep.o matrix_functions.o nthLargestPM.o preprocess.o psi_fns.o qnorm_probeset.o read_rmaexpress.o rlm_PLM.o rlm_threestep.o rmaPLM_pseudo.o rma_PLM.o rma_common.o scaling.o threestep.o threestep_PLM.o threestep_common.o threestep_summary.o threestep_summary_methods.o transfns.o -L/home/biocbuild/bbs-3.11-bioc/R/lib -lRlapack -L/home/biocbuild/bbs-3.11-bioc/R/lib -lRblas -lgfortran -lm -lquadmath -lz -L/home/biocbuild/bbs-3.11-bioc/R/lib -lR
installing to /home/biocbuild/bbs-3.11-bioc/R/library/00LOCK-affyPLM/00new/affyPLM/libs
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** checking absolute paths in shared objects and dynamic libraries
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (affyPLM)
affyPLM.Rcheck/tests/C_code_tests.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ####
> #### This code is messy, possibly incomplete and only for
> #### the use of developers.
> ####
> ####
>
> test.c.code <- FALSE
> test.PLM.modelmatrix <- FALSE
> test.rlm <- FALSE
>
> if (test.c.code){
+
+ library(affyPLM)
+ narrays <- 10
+ nprobes <- 11
+ nprobetypes <- 2
+ ncols <- 10
+
+ MMs <- rnorm(narrays*nprobes*nprobetypes)
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+
+ #test making intercept column
+ matrix(.C("R_PLM_matrix_intercept",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),0)[[1]],ncol=10)
+
+ #test making an MM covariate column
+ matrix(.C("R_PLM_matrix_MM",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.double(MMs))[[1]],ncol=10)
+
+ # sample effect aka chip effect, aka expression values
+ matrix(.C("R_PLM_matrix_sample_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0))[[1]],ncol=10)
+ matrix(.C("R_PLM_matrix_sample_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1))[[1]],ncol=10)
+ matrix(.C("R_PLM_matrix_sample_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1))[[1]],ncol=10)
+
+
+
+ #probe-type parameter overall
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(0),integer(narrays),as.integer(1))[[1]],ncol=10)
+
+ #probe-type parameter within sample
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(1),integer(narrays),as.integer(1))[[1]],ncol=10)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(1),integer(narrays),as.integer(1))[[1]],ncol=10)
+ ncols <- 20
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(1),integer(narrays),as.integer(0))[[1]],ncol=20)
+
+
+ #probe-type-parameter within a chip-level factor (eg treatment, or genotype variable)
+ trt.cov <- rep(0:1,5)
+ ncols <- 10
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=10)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=10)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=10)
+
+ trt.cov <- rep(0:4,2)
+ matrix(.C("R_PLM_matrix_probe_type_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(2),as.integer(trt.cov),as.integer(4))[[1]],ncol=10)
+
+
+
+
+ #probe effects - overall
+ ncols <- 11
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+
+
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(0),as.integer(trt.cov),as.integer(4))[[1]],ncol=11)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(0),as.integer(trt.cov),as.integer(4))[[1]],ncol=11)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(0),as.integer(trt.cov),as.integer(4))[[1]],ncol=11)
+
+
+
+ #probe effects within treatment or genotype factor
+ trt.cov <- rep(0:1,5)
+ ncols <- 22
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(2),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+
+
+ #probe effects within probetype
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(3),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(3),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(3),as.integer(trt.cov),as.integer(1))[[1]],ncol=22)
+
+
+ #probe effects within probetype within treatment or genotype factor variable
+ trt.cov <- rep(0:1,5)
+ ncols <- 44
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+
+ nprobetypes <- 1
+ trt.cov <- rep(0:1,5)
+ ncols <- 44
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+ matrix(.C("R_PLM_matrix_probe_effect",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(-1),as.integer(4),as.integer(trt.cov),as.integer(1))[[1]],ncol=44)
+
+
+ # copy across chip level variables into model matrix
+ nprobetypes <- 1
+ trt.cov <- rep(0:1,5)
+ ncols <- 10
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ trt.variables <- rnorm(10)
+
+ matrix(.C("R_PLM_matrix_chiplevel",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.double(trt.variables),as.integer(1))[[1]],ncol=10)
+
+
+ ###
+ ### Build a few design matrices and compare with R model.matrix
+ ###
+
+
+ for (nprobetypes in 1:2){
+ for (narrays in 2:15){
+ for (nprobes in 2:20){
+ for (constraint.type in c("contr.sum","contr.treatment")){
+ if (constraint.type == "contr.sum"){
+ ct.type <- -1
+ } else {
+ ct.type <- 1
+ }
+
+
+ ncols <- nprobes -1 + narrays
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(1),as.integer(0),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ probe.effect <- factor(rep(1:nprobes,narrays*nprobetypes))
+ sample.effect <- factor(rep(rep(c(1:narrays),rep(nprobes,narrays)),nprobetypes))
+ if (nprobetypes == 2){
+ probe.type.effect <- factor(rep(1:2,c(narrays*nprobes,narrays*nprobes)))
+ } else {
+ probe.type.effect <- factor(rep(1,narrays*nprobes))
+ }
+
+ if (any(X!=model.matrix(~ C(sample.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes, " ", nprobetypes)
+ }
+
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(0),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ -1 + C(sample.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- nprobes
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(0),as.integer(0),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+
+ ncols <- nprobes
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(0),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~-1+ C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ }
+ }
+ }
+ }
+
+ ###
+ ### Build a few more design matrices and compare with R model.matrix
+ ###
+
+
+ for (narrays in 2:15){
+ for (nprobes in 2:20){
+ for (constraint.type in c("contr.sum","contr.treatment")){
+ probe.effect <- factor(rep(1:nprobes,narrays*nprobetypes))
+ sample.effect <- factor(rep(rep(c(1:narrays),rep(nprobes,narrays)),nprobetypes))
+ if (constraint.type == "contr.sum"){
+ ct.type <- -1
+ } else {
+ ct.type <- 1
+ }
+
+
+ if (nprobetypes == 2){
+ probe.type.effect <- factor(rep(1:2,c(narrays*nprobes,narrays*nprobes)))
+ } else {
+ probe.type.effect <- factor(rep(1,narrays*nprobes))
+ }
+ ncols <- nprobetypes + nprobes -1
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(1),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~-1+ C(probe.type.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- nprobetypes + nprobes -1
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(0),as.integer(1),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ C(probe.type.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- narrays + nprobetypes + nprobes -2
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(1),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ -1 + C(sample.effect,constraint.type) + C(probe.type.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- narrays + nprobetypes + nprobes -2
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(1),as.integer(1),as.integer(1),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ + C(sample.effect,constraint.type) + C(probe.type.effect,constraint.type) + C(probe.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- narrays + nprobetypes -1
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(1),as.integer(1),as.integer(0),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ C(sample.effect,constraint.type) + C(probe.type.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- narrays + nprobetypes -1
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(1),as.integer(1),as.integer(0),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ -1 + C(sample.effect,constraint.type) + C(probe.type.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- nprobetypes
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(1),as.integer(0),as.integer(1),as.integer(0),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~ C(probe.type.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ ncols <- nprobetypes
+ X <- matrix(0,narrays*nprobes*nprobetypes,ncols)
+ X <- matrix(.C("R_PLM_Matrix_constructtest",as.double(as.vector(X)), as.integer(narrays),as.integer(nprobes),as.integer(nprobetypes),as.integer(0),as.integer(0),as.integer(1),as.integer(0),as.integer(ct.type))[[1]],ncol=ncols)
+
+ if (any(X!=model.matrix(~-1 + C(probe.type.effect,constraint.type)))){
+ stop("Model matrix function problem ",narrays," ", nprobes," ", nprobetypes)
+ }
+
+ }
+
+
+ }
+ }
+
+
+ narrays <- 2
+ nprobes <- 7
+ nprobetypes <- 2
+
+ probe.effect <- factor(rep(1:nprobes,narrays*nprobetypes))
+ sample.effect <- factor(rep(rep(c(1:narrays),rep(nprobes,narrays)),nprobetypes))
+ if (constraint.type == "contr.sum"){
+ ct.type <- -1
+ } else {
+ ct.type <- 1
+ }
+
+
+ if (nprobetypes == 2){
+ probe.type.effect <- factor(rep(1:2,c(narrays*nprobes,narrays*nprobes)))
+ } else {
+ probe.type.effect <- factor(rep(1,narrays*nprobes))
+ }
+
+
+ model.matrix(~-1 +probe.effect/probe.type.effect)
+
+
+ library(affyPLM)
+ output <- verify.output.param(list(weights = FALSE, residuals = FALSE, varcov = "none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ library(affydata)
+ data(Dilution)
+
+ # fit a PM ~ samples model
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0), probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ sample.effect <- rep(1:4,c(16,16,16,16))
+ probe.effect <- rep(1:16,4)
+
+ library(MASS)
+ fit <- rlm(as.vector(log2(pm(Dilution)[1:16,])) ~ -1 + factor(sample.effect))
+
+ if (any(Fitresults[[1]][1,] != coef(fit))){
+ stop("Problem in model fitting procedure")
+ }
+
+ sample.effect <- rep(1:4,c(20,20,20,20))
+ probe.effect <- rep(1:20,4)
+ fit <- rlm(as.vector(log2(pm(Dilution)[201781:201800,])) ~ -1 + factor(sample.effect))
+ if (any(Fitresults[[1]][12625,] != coef(fit))){
+ stop("Problem in model fitting procedure")
+ }
+
+
+ # fit a samples + probes model
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0), probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+ sample.effect <- rep(1:4,c(16,16,16,16))
+ probe.effect <- rep(1:16,4)
+
+ fit <- rlm(as.vector(log2(pm(Dilution)[1:16,])) ~ -1 + factor(sample.effect)+C(factor(probe.effect),"contr.sum"))
+
+ if (any(abs(Fitresults[[1]][1,] -coef(fit)[1:4]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[2]][[1]]) - coef(fit)[5:19]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+
+
+ sample.effect <- rep(1:4,c(20,20,20,20))
+ probe.effect <- rep(1:20,4)
+ fit <- rlm(as.vector(log2(pm(Dilution)[201781:201800,])) ~ -1 + factor(sample.effect)+C(factor(probe.effect),"contr.sum"))
+
+ if (any(abs(Fitresults[[1]][12625,] -coef(fit)[1:4])> 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[2]][[12625]])- coef(fit)[5:23])>1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ # fit an MM ~ samples model
+ R.model <- list(mmorpm.covariate=0,response.variable=-1,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0), probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ sample.effect <- rep(1:4,c(16,16,16,16))
+ probe.effect <- rep(1:16,4)
+
+ library(MASS)
+ fit <- rlm(as.vector(log2(mm(Dilution)[1:16,])) ~ -1 + factor(sample.effect))
+
+ if (any(abs(Fitresults[[1]][1,] - coef(fit)) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ sample.effect <- rep(1:4,c(20,20,20,20))
+ probe.effect <- rep(1:20,4)
+ fit <- rlm(as.vector(log2(mm(Dilution)[201781:201800,])) ~ -1 + factor(sample.effect))
+ if (any(abs(Fitresults[[1]][12625,] - coef(fit))>1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ # fit a MM ~ samples + probes model
+ R.model <- list(mmorpm.covariate=0,response.variable=-1,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0), probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+ sample.effect <- rep(1:4,c(16,16,16,16))
+ probe.effect <- rep(1:16,4)
+
+ fit <- rlm(as.vector(log2(mm(Dilution)[1:16,])) ~ -1 + factor(sample.effect)+C(factor(probe.effect),"contr.sum"))
+
+ if (any(abs(Fitresults[[1]][1,]-coef(fit)[1:4]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[2]][[1]])- coef(fit)[5:19])>1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+
+
+ sample.effect <- rep(1:4,c(20,20,20,20))
+ probe.effect <- rep(1:20,4)
+ fit <- rlm(as.vector(log2(mm(Dilution)[201781:201800,])) ~ -1 + factor(sample.effect)+C(factor(probe.effect),"contr.sum"))
+
+ if (any(abs(Fitresults[[1]][12625,]- coef(fit)[1:4])>1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[2]][[12625]])-coef(fit)[5:23])>1e14)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+ # a treatment model
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- model.matrix(~ -1 + as.factor(treatment.effect))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=-1,which.parameter.types=as.integer(c(0,1,0,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =covariates, probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ treatment.effect <- rep(c(1,1,2,2),c(16,16,16,16))
+ fit <- rlm(as.vector(log2(mm(Dilution)[1:16,])) ~ -1 + factor(treatment.effect))
+
+ if (any(abs(Fitresults[[1]][1,]-coef(fit)[1:2]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ output <- verify.output.param(list(weights = FALSE, residuals = FALSE, varcov = "none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+
+ # a treatment + probes model with contr.treatment constraint
+ R.model <- list(mmorpm.covariate=0,response.variable=-1,which.parameter.types=as.integer(c(0,1,0,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =covariates, probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ treatment.effect <- rep(c(1,1,2,2),c(20,20,20,20))
+ probe.effect <- rep(1:20,4)
+ fit <- rlm(as.vector(log2(mm(Dilution)[201761:201780,])) ~ -1 + factor(treatment.effect)+C(factor(probe.effect),"contr.treatment"))
+
+ if (any(abs(Fitresults[[1]][12624,]-coef(fit)[1:2]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[2]][[12624]])-coef(fit)[3:21])>1e14)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+
+
+ # MM + samples + probes
+ R.model <- list(mmorpm.covariate=1,response.variable=1,which.parameter.types=as.integer(c(0,0,0,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0), probe.type.levels=list(),probe.trt.levels=list())
+
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ sample.effect <- rep(1:4,c(16,16,16,16))
+ probe.effect <- rep(1:16,4)
+
+ fit <- rlm(as.vector(log2(pm(Dilution)[1:16,])) ~ -1 + as.vector(log2(mm(Dilution)[1:16,])))
+
+
+ if (any(abs(as.vector(Fitresults[[6]][[1]]) - coef(fit)) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ R.model <- list(mmorpm.covariate=1,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ fit <- rlm(as.vector(log2(pm(Dilution)[1:16,])) ~ -1 + as.vector(log2(mm(Dilution)[1:16,]))+ as.factor(sample.effect))
+ if (any(abs(as.vector(Fitresults[[6]][[1]]) - coef(fit)[1]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[1]][1,]) - coef(fit)[2:5]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ R.model <- list(mmorpm.covariate=1,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ fit <- rlm(as.vector(log2(pm(Dilution)[1:16,])) ~ -1 + as.vector(log2(mm(Dilution)[1:16,]))+ as.factor(sample.effect) + C(as.factor(probe.effect),"contr.sum"))
+ if (any(abs(as.vector(Fitresults[[6]][[1]]) - coef(fit)[1]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[1]][1,]) - coef(fit)[2:5]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+ ## PM and MM are response
+
+
+ sample.effect <- rep(1:4,c(32,32,32,32))
+ probe.effect <- rep(1:16,8)
+
+
+
+ # PMMM ~ -1 + samples
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ fit <- rlm(as.vector(rbind(log2(pm(Dilution)[1:16,]),log2(mm(Dilution)[1:16,]))) ~ -1 + as.factor(sample.effect))
+ if (any(abs(as.vector(Fitresults[[1]][1,]) - coef(fit)) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ # PMMM ~ -1 + samples +PROBES
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ fit <- rlm(as.vector(rbind(log2(pm(Dilution)[1:16,]),log2(mm(Dilution)[1:16,]))) ~ -1 + as.factor(sample.effect)+C(as.factor(probe.effect),"contr.sum") )
+ if (any(abs(as.vector(Fitresults[[1]][1,]) - coef(fit)[1:4]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ # a probe.type effect
+ probe.type.effect <- rep(rep(1:2,c(16,16)),4)
+
+ # PMMM ~ -1 + samples + probe.type + PROBES
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,1,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,-1,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ fit <- rlm(as.vector(rbind(log2(pm(Dilution)[1:16,]),log2(mm(Dilution)[1:16,]))) ~ -1 + as.factor(sample.effect)+ C(as.factor(probe.type.effect),"contr.sum")+ C(as.factor(probe.effect),"contr.sum") )
+
+ if (any(abs(as.vector(Fitresults[[1]][1,]) - coef(fit)[1:4]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+ if (any(abs(as.vector(Fitresults[[6]][1,]) - coef(fit)[5]) > 1e-13)){
+ stop("Problem in model fitting procedure")
+ }
+
+
+
+ #### store weights PM ~ -1 + samples
+
+ output <- verify.output.param(list(weights = TRUE, residuals = TRUE, varcov = "none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ #### store weights PMMM ~ -1 + samples
+
+ output <- verify.output.param(list(weights = TRUE, residuals = TRUE, varcov = "none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ #### store weights PMMM ~ -1 + samples + probe.type + probes
+
+ output <- verify.output.param(list(weights = TRUE, residuals = TRUE, varcov ="none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,1,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,1,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PM ~ -1 + treatment + probes in treatment
+ output <- verify.output.param(list(weights = TRUE, residuals = TRUE, varcov = "none", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- model.matrix(~ -1 + as.factor(treatment.effect))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,1,0,0,1)),strata=as.integer(c(0,0,0,0,2)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =covariates,probe.type.levels=list(),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1 + treatment + probes in treatment
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- model.matrix(~ -1 + as.factor(treatment.effect))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,1,0,0,1)),strata=as.integer(c(0,0,0,0,2)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =covariates,probe.type.levels=list(),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1 + treatment + probe.effect in treatment
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- model.matrix(~ -1 + as.factor(treatment.effect))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,1,0,1,0)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,-1,0)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =covariates,probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,0)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes.type + probes
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=NULL,max.probe.type.trt.factor=0,probe.trt.factor=NULL,max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list(),probe.trt.levels=list())
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes.type + probes with both within treatment factor
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,1)),strata=as.integer(c(0,0,0,2,2)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ ## PMMM ~ -1+ probes.type + probes with both within treatment factor and probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,1)),strata=as.integer(c(0,0,0,2,4)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ ## PMMM ~ -1+ probes.type + probes probe.types within treatment factor and probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,1)),strata=as.integer(c(0,0,0,2,3)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes.type + probes probe.types within samples and probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,1,1)),strata=as.integer(c(0,0,0,1,3)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ ## PMMM ~ intercept + probes.type + probes probe.types within samples and probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,0,0,1,1)),strata=as.integer(c(0,0,0,1,3)),constraints=as.integer(c(0,0,0,-1,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ intercept + probes probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,0,0,0,1)),strata=as.integer(c(0,0,0,0,3)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## PMMM ~ -1+ probes probes also within probe.type
+ output <- list(weights = TRUE, residuals = TRUE, varcov = c("none","chiplevel", "all"), resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,0,0,1)),strata=as.integer(c(0,0,0,0,3)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ # now play with varcov output
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ # now play with varcov output and treatment
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- model.matrix(~ -1 + as.factor(treatment.effect))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,1,0,0,0)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =covariates,probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ # now play with varcov output and an intercept
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,0,1,0,1)),strata=as.integer(c(0,0,0,0,0)),constraints=as.integer(c(0,0,-1,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=1,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ # now play with varcov output and treatment and intercept
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- matrix(model.matrix(~ as.factor(treatment.effect))[,2])
+ colnames(covariates) <- "trt_2"
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,1,0,0,0)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =covariates,probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ # now play with varcov all option output and treatment and intercept
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ treatment.effect <- c(1,1,2,2)
+
+ covariates <- matrix(model.matrix(~ as.factor(treatment.effect))[,2])
+ colnames(covariates) <- "trt_2"
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,1,0,0,0)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,0)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =covariates,probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ # now play with varcov all option output and samples and intercept
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(1,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,-1,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ # now play with varcov all option output and samples and intercept, MM covarite
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=1,response.variable=1,which.parameter.types=as.integer(c(1,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,-1,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+
+
+
+
+ ## now play with varcov all option output and samples and intercept, MM covariate and input chip weights
+ output <- verify.output.param(list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE))
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=c(1,1,0.5,0.5),weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ## now play with varcov all option output and samples and intercept, MM covariate and input chip weights
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=runif(201800))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=1,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=c(rep(c(1,0.5),c(201800,201800))))
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="cuberoot", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log10", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+
+ R.model <- list(mmorpm.covariate=0,response.variable=0,which.parameter.types=as.integer(c(0,0,1,0,1)),strata=as.integer(c(0,0,0,2,0)),constraints=as.integer(c(0,0,0,0,-1)),probe.type.trt.factor=as.integer(c(0,0,1,1)),max.probe.type.trt.factor=1,probe.trt.factor=as.integer(c(0,0,1,1)),max.probe.trt.factor=0,chipcovariates =matrix(0,0,0),probe.type.levels=list("blah"=c("A","B")),probe.trt.levels=list("blah"=c("A","B")))
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ }
>
>
>
> if (test.PLM.modelmatrix){
+
+ library(affyPLM);data(Dilution)
+
+ #PLM.designmatrix3(Dilution)
+
+ #PLM.designmatrix3(Dilution,model=MM ~ PM -1 + samples +probe.type:probes)
+
+ #PLM.designmatrix3(Dilution,model=MM ~ PM -1 + samples:probe.type + liver:probe.type:probes + liver:samples)
+ #PLM.designmatrix3(Dilution,model=MM ~ PM + samples:probe.type + liver:probe.type:probes + liver + samples)
+
+
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ #blah <- c(1,5,5,1)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ probes + blah,constraint.type=c(probes="contr.sum"))
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 + blah:probe.type)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 +probes:probe.type)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 +probes:blah)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 +probes:probe.type:blah)
+ #output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ # R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ samples,constraint.type=c(samples="contr.sum"))
+ # R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ blah)
+ # R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 + samples)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ probes + blah)
+ #R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ probes + blah)
+ #Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ library(affyPLM);data(Dilution)
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ probe.type + probe.type:probes + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ library(affyPLM);data(Dilution)
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ samples:probe.type + probe.type:probes + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="chiplevel", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ blah <- c(1,2,2)
+ R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ blah:probe.type + probe.type:probes + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ blah <- c(1,2,2)
+ R.model <- PLM.designmatrix3(Dilution,model=PM ~ -1 + probes + MM + blah)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+ output <- list(weights = TRUE, residuals = TRUE, varcov ="all", resid.SE = TRUE)
+ modelparam <- list(trans.fn="log2", se.type = 4, psi.type = 0, psi.k =1.345,max.its = 20, init.method = "ls",weights.chip=NULL,weights.probe=NULL)
+ blah <- c(1,2,2)
+ R.model <- PLM.designmatrix3(Dilution,model=PM ~ -1 + probes + MM + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+
+
+ #test some of the verification functions
+
+
+ output <- verify.output.param()
+ modelparam <- verify.model.param(Dilution,PM ~ -1 + probes + MM + samples)
+ R.model <- PLM.designmatrix3(Dilution,model=PM ~ -1 + probes + MM + samples)
+
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+ ##verify.model.param(Dilution,PM ~ -1 + probes + MM + samples,model.param=list(weights.probe=rep(1,10)))
+
+ modelparam <- verify.model.param(Dilution,PMMM ~ -1 + probes + samples,model.param=list(weights.chip=c(1,2,3),weights.probe=rep(1,2400*2)))
+ R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 + probes + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ modelparam <- verify.model.param(Dilution,PM ~ -1 + probes + samples,model.param=list())
+ R.model <- PLM.designmatrix3(Dilution,model=PM ~ -1 + probes + samples)
+ Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+ ## probes <- rep(1:16,3)
+ ## chips <- rep(1:3,c(16,16,16))
+
+ ## library(MASS)
+
+ ##fit <- rlm(log2(as.vector(pm(Dilution,"HG2188-HT2258_at"))) ~ -1 + as.factor(chips) + C(as.factor(probes),"contr.sum"))
+
+
+ #test creating a PLMset based on the output from rlm_PLMset
+
+ ### x <- new("PLMset")
+ ### x@chip.coefs=Fitresults[[1]]
+ ### x@probe.coefs= Fitresults[[2]]
+ ### x@weights=Fitresults[[3]]
+ ### x@se.chip.coefs=Fitresults[[4]]
+ ### x@se.probe.coefs=Fitresults[[5]]
+ ### x@exprs=Fitresults[[6]]
+ ### x@se.exprs=Fitresults[[7]]
+ ### x@residuals=Fitresults[[8]]
+ ### x@residualSE=Fitresults[[9]]
+ ### x@varcov = Fitresults[[10]]
+ ### x@cdfName = Dilution@cdfName
+ ### x@phenoData = Dilution@phenoData
+ ### x@annotation = Dilution@annotation
+ ### x@description = Dilution@description
+ ### x@notes = Dilution@notes
+ ### x@nrow= Dilution@nrow
+ ### x@ncol= Dilution@ncol
+ ### x@model.description = c(x@model.description, list(R.model=R.model))
+ ### image(x)
+
+
+
+
+ ### data(Dilution)
+ ### output <- verify.output.param()
+ ### modelparam <- verify.model.param(Dilution,PMMM ~ -1 + probe.type:probes + samples + samples:probe.type,model.param=list())
+ ### R.model <- PLM.designmatrix3(Dilution,model=PMMM ~ -1 + probe.type:probes + samples+ samples:probe.type)
+ ### Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+ ### output <- verify.output.param()
+ ### modelparam <- verify.model.param(Dilution,MM ~ -1 + probes + samples,model.param=list())
+ ### R.model <- PLM.designmatrix3(Dilution,model=MM ~ -1 + probes + samples)
+ ### Fitresults <- .Call("rlm_PLMset",pm(Dilution),mm(Dilution),probeNames(Dilution),length(geneNames(Dilution)),R.model,output,modelparam)
+
+
+
+ ### x <- new("PLMset")
+ ### x@chip.coefs=Fitresults[[1]]
+ ### x@probe.coefs= Fitresults[[2]]
+ ### x@weights=Fitresults[[3]]
+ ### x@se.chip.coefs=Fitresults[[4]]
+ ### x@se.probe.coefs=Fitresults[[5]]
+ ### x@exprs=Fitresults[[6]]
+ ### x@se.exprs=Fitresults[[7]]
+ ### x@residuals=Fitresults[[8]]
+ ### x@residualSE=Fitresults[[9]]
+ ### x@varcov = Fitresults[[10]]
+ ### x@cdfName = Dilution@cdfName
+ ### x@phenoData = Dilution@phenoData
+ ### x@annotation = Dilution@annotation
+ ### x@description = Dilution@description
+ ### x@notes = Dilution@notes
+ ### x@nrow= Dilution@nrow
+ ### x@ncol= Dilution@ncol
+ ### x@model.description = c(x@model.description, list(R.model=R.model))
+ ### image(x)
+ ### image(x,type="pos.resids")
+ ### image(x,type="neg.resids")
+ ### image(x,type="sign.resids")
+
+ ### resid(x,"1091_at")
+
+
+
+ ### weights(x,c("1091_at","1092_at"))
+
+
+ ### image(x,type="resids",standardize=TRUE)
+
+
+
+
+
+
+ }
>
>
>
>
>
> if (test.rlm){
+
+
+ library(affyPLM);data(Dilution)
+
+ y <- as.vector(log2(pm(Dilution)[1:16,]))
+
+ w <- runif(64)
+
+ probes <- rep(1:16,4)
+ samples <- rep(1:4,c(16,16,16,16))
+
+ x <- model.matrix( ~ -1 + as.factor(samples) + C(as.factor(probes),"contr.sum"))
+ x <- as.vector(x)
+
+ cols <- 19
+ rows <- 64
+
+
+ # rlm_wfit_R(double *x, double *y, double *w, int *rows, int *cols, double *out_beta, double *out_resids, double *out_weights)
+
+ fit1 <-.C("rlm_wfit_R",as.double(x),as.double(y),as.double(w),as.integer(rows),as.integer(cols),double(cols),double(rows),double(rows))
+
+
+ library(MASS)
+
+ fit2 <- rlm(y ~ -1 + as.factor(samples) + C(as.factor(probes),"contr.sum"),weights=w,wt.method="case")
+
+ if (any(abs(coef(fit2) - fit1[[6]]) > 10e-14)){
+ stop("Weighted RLM did not work")
+ }
+
+
+
+
+
+ y <- as.vector(log2(pm(Dilution,"1001_at")))
+ x <- as.vector(log2(mm(Dilution,"1001_at")))
+
+ rlm(y ~ -1 + x + as.factor(samples) + C(as.factor(probes),"contr.sum"))
+
+
+
+
+
+
+
+
+
+
+
+ }
>
>
>
>
> proc.time()
user system elapsed
3.713 0.055 3.754
affyPLM.Rcheck/tests/PLM_tests.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> do.all.tests <- FALSE
> if (do.all.tests){
+
+ # this file tests fitPLM and the PLMset object
+
+ library(affyPLM)
+
+ library(affydata)
+ data(Dilution)
+
+
+ Pset <- fitPLM(Dilution)
+
+ #check accessors for parameters and se
+
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+ coefs.probe(Pset)[1:5]
+ se.probe(Pset)[1:5]
+ coefs.const(Pset)
+ se.const(Pset)
+
+ #accessors for weights and residuals
+
+ weights(Pset)[[1]][1:5,]
+ resid(Pset)[[1]][1:5,]
+
+
+ #test varcov
+
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,output.param=list(varcov="chiplevel"))
+ varcov(Pset)[1:3]
+
+
+ #test each of the possible weight functions
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Huber"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="fair"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Cauchy"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Geman-McClure"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Welsch"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Tukey"))
+ Pset <- fitPLM(Dilution,background=FALSE,normalize=FALSE,model.param=list(psi.type="Andrews"))
+
+ # a larger example to do some testing of the graphical functions
+
+ data(Dilution)
+
+ Pset <- fitPLM(Dilution)
+
+ #testing the image capabilities
+
+ image(Pset,which=2)
+ image(Pset,which=2,type="resids")
+ image(Pset,which=2,type="pos.resids")
+ image(Pset,which=2,type="neg.resids")
+ image(Pset,which=2,type="resids",use.log=FALSE,add.legend=TRUE)
+
+ boxplot(Pset)
+ Mbox(Pset)
+
+
+ #test some non-default models functions
+ # no preprocessing for speed
+
+ Pset <- fitPLM(Dilution, PM ~ -1 + probes + liver,background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+
+
+ Pset <- fitPLM(Dilution, PM ~ -1 + probes + liver + scanner,background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+
+ #checking the constraints
+ Pset <- fitPLM(Dilution, PM ~ -1 + probes + liver + scanner,constraint.type=c(default="contr.sum"),background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+
+ Pset <- fitPLM(Dilution, PM ~ -1 + liver + scanner,constraint.type=c(default="contr.sum"),background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+ coefs.probe(Pset) # should be empty
+
+ Pset <- fitPLM(Dilution, PM ~ -1 + liver + scanner,constraint.type=c(probes="contr.treatment"),background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+ coefs.probe(Pset) # should be empty
+
+
+ Pset <- fitPLM(Dilution, PM ~ -1 + probes + liver + scanner,constraint.type=c(probes="contr.treatment"),background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+ coefs.probe(Pset)[1:16]
+
+
+ scanner2 <- c(1,2,1,2)
+ Pset <- fitPLM(Dilution, PM ~ -1 + probes + liver + scanner2,constraint.type=c(probes="contr.sum"),background=FALSE,normalize=FALSE)
+ coefs(Pset)[1:5,]
+ se(Pset)[1:5,]
+ coefs.probe(Pset)[1:16]
+
+ #
+ #Pset <- fitPLM(Dilution,model=PM~-1+probes+scanner,normalize=FALSE,background=FALSE,model.param=list(se.type=3))
+ #se(Pset)[1:10,]
+
+ #check that fitPLM rlm agrees with threestep rlm and threestepPLM rlm
+
+
+ Pset <- fitPLM(Dilution)
+ eset <- threestep(Dilution,summary.method="rlm")
+ Pset2 <- threestepPLM(Dilution,summary.method="rlm")
+
+ if (any(abs(coefs(Pset) - exprs(eset)) > 1e-14)){
+ stop("no agreement between fitPLM and threestep")
+ }
+
+ if (any(abs(coefs(Pset) - coefs(Pset2)) > 1e-14)){
+ stop("no agreement between fitPLM and threestep")
+ }
+ }
>
> proc.time()
user system elapsed
0.306 0.014 0.305
affyPLM.Rcheck/tests/preprocess_tests.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> #test the preprocessing functionality
>
> library(affyPLM)
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':
Filter, Find, Map, Position, Reduce, anyDuplicated, append,
as.data.frame, basename, cbind, colnames, dirname, do.call,
duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
tapply, union, unique, unsplit, which, which.max, which.min
Loading required package: affy
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: gcrma
Loading required package: preprocessCore
> library(affydata)
Package LibPath Item
[1,] "affydata" "/home/biocbuild/bbs-3.11-bioc/R/library" "Dilution"
Title
[1,] "AffyBatch instance Dilution"
> data(Dilution)
>
>
> ### NO LONGER SUPPORTED eset <- threestep(Dilution,background.method="RMA.1")
> eset <- threestep(Dilution,background.method="RMA.2")
Warning messages:
1: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu95av2cdf'
2: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu95av2cdf'
> eset <- threestep(Dilution,background.method="IdealMM")
> eset <- threestep(Dilution,background.method="MAS")
> eset <- threestep(Dilution,background.method="MASIM")
> eset <- threestep(Dilution,background.method="LESN2")
> eset <- threestep(Dilution,background.method="LESN1")
> eset <- threestep(Dilution,background.method="LESN0")
>
> eset <- threestep(Dilution,normalize.method="quantile",background=FALSE)
> eset <- threestep(Dilution,normalize.method="quantile.probeset",background=FALSE)
> eset <- threestep(Dilution,normalize.method="scaling",background=FALSE)
>
>
>
> proc.time()
user system elapsed
22.929 0.386 23.316
affyPLM.Rcheck/tests/threestepPLM_tests.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> if (.Platform$OS.type != "windows"){
+ library(affyPLM)
+
+ # test threestep and threestepPLM to see if they agree
+
+
+ check.coefs <- function(Pset,Pset2){
+ if (any(abs(coefs(Pset) - exprs(Pset2)) > 1e-14)){
+ stop("No agreement between threestepPLM and threestep in coefs")
+ }
+ }
+
+ check.resids <- function(Pset,Pset2){
+ if (any(resid(Pset) != resid(Pset2))){
+ stop("No agreement between threestepPLM and rmaPLM/threestep in residuals")
+ }
+ }
+
+
+ library(affydata)
+ data(Dilution)
+
+ Pset <- threestepPLM(Dilution)
+ Pset2 <- threestep(Dilution)
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="tukey.biweight")
+ Pset2 <- threestep(Dilution,summary.method="tukey.biweight")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="average.log")
+ Pset2 <- threestep(Dilution,summary.method="average.log")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="rlm")
+ Pset2 <- threestep(Dilution,summary.method="rlm")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="log.average")
+ Pset2 <- threestep(Dilution,summary.method="log.average")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="log.median")
+ Pset2 <- threestep(Dilution,summary.method="log.median")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="median.log")
+ Pset2 <- threestep(Dilution,summary.method="median.log")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="log.2nd.largest")
+ Pset2 <- threestep(Dilution,summary.method="log.2nd.largest")
+ check.coefs(Pset,Pset2)
+
+ Pset <- threestepPLM(Dilution,summary.method="lm")
+ Pset2 <- threestep(Dilution,summary.method="lm")
+ check.coefs(Pset,Pset2)
+
+ #check if threestepPLM agrees with rmaPLM
+ Pset <- threestepPLM(Dilution)
+ Pset2 <- rmaPLM(Dilution)
+
+ if (any(coefs(Pset) != coefs(Pset2))){
+ stop("No agreement between threestepPLM and rmaPLM in coefs")
+ }
+
+
+ if (any(resid(Pset)[[1]] != resid(Pset2)[[1]])){
+ stop("No agreement between threestepPLM and rmaPLM in residuals")
+ }
+ }
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':
Filter, Find, Map, Position, Reduce, anyDuplicated, append,
as.data.frame, basename, cbind, colnames, dirname, do.call,
duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
tapply, union, unique, unsplit, which, which.max, which.min
Loading required package: affy
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: gcrma
Loading required package: preprocessCore
Package LibPath Item
[1,] "affydata" "/home/biocbuild/bbs-3.11-bioc/R/library" "Dilution"
Title
[1,] "AffyBatch instance Dilution"
Warning messages:
1: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu95av2cdf'
2: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu95av2cdf'
>
> proc.time()
user system elapsed
59.658 0.225 59.868
affyPLM.Rcheck/affyPLM-Ex.timings
| name | user | system | elapsed | |
| PLMset2exprSet | 6.726 | 0.119 | 6.847 | |
| bg.correct.LESN | 1.207 | 0.023 | 1.232 | |
| fitPLM | 12.353 | 0.155 | 12.508 | |
| normalize.exprSet | 1.110 | 0.019 | 1.131 | |
| normalize.scaling | 1.533 | 0.031 | 1.566 | |
| preprocess | 2.147 | 0.008 | 2.155 | |
| rmaPLM | 0.506 | 0.007 | 0.514 | |
| threestep | 18.620 | 0.007 | 18.631 | |
| threestepPLM | 0.314 | 0.004 | 0.317 | |