[Genabel-commits] r1967 - in pkg/MixABEL: R tests tests/broken
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 10 09:44:20 CEST 2015
Author: lckarssen
Date: 2015-05-10 09:44:20 +0200 (Sun, 10 May 2015)
New Revision: 1967
Added:
pkg/MixABEL/tests/broken/
pkg/MixABEL/tests/broken/rwth_caller.R
pkg/MixABEL/tests/broken/test_rwth.R
Removed:
pkg/MixABEL/R/rwth_caller.R
pkg/MixABEL/tests/test_rwth.R
Log:
Moved broken MixABEL unittests to a separate directory.
Deleted: pkg/MixABEL/R/rwth_caller.R
===================================================================
--- pkg/MixABEL/R/rwth_caller.R 2015-05-08 13:00:25 UTC (rev 1966)
+++ pkg/MixABEL/R/rwth_caller.R 2015-05-10 07:44:20 UTC (rev 1967)
@@ -1,60 +0,0 @@
-"rwth_caller" <- function(
- printall,
-# parameters and data to compute M
- Phi, H2, Sigma2,
-# vector of length length_y
- y,
-# matrix of height length_y and width widthXL
- XL,
-# matrix of height length_y and width nXRs*widthXRs
-# it contains nXRs matrices XR
- XRs,
-# width of single XR
- WidthXR
-)
-{
-# figure out data dims and some sanity checks
- Length_y = length(y);
- NXRs = dim(XRs)[2]/WidthXR;
- if ((dim(XRs)[2] %% WidthXR) != 0) stop("(dim(XRs)[2] %% widthXR) != 0");
- if (dim(XRs)[1] != length(y)) stop("dim(XRs)[1] != length(y)");
- if (is.null(XL)) {
- XL <- matrix(1,ncol=1,nrow=length(y));
- } else {
- tmp <- matrix(1,ncol=1,nrow=length(y));
- XL <- cbind(tmp,XL)
- }
- if (dim(XL)[1] != length(y)) stop("dim(XL)[1] != length(y)");
- WidthXL = dim(XL)[2];
-
- ncol_beta <- WidthXL + WidthXR;
- nrow_out <- NXRs;
- ncol_vcbeta <- ncol_beta*(ncol_beta+1)/2
-
- out <- .C("rwth_example",
- as.integer(printall),
-# dimensions of the data, explained later
- as.integer(Length_y), as.integer(NXRs),
- as.integer(WidthXR), as.integer(WidthXL),
- # parameters and data to compute M
- as.double(Phi), as.double(H2), as.double(Sigma2),
- # vector of length length_y
- as.double(y),
- # matrix of height length_y and width widthXL
- as.double(XL),
- # matrix of height length_y and width nXRs*widthXRs
- # it contains nXRs matrices XR
- as.double(XRs),
- # this is output
- # 'beta's length is (widthXL + widthXRs)
- # for time-being, just fill it with mean of
- # columns of X
- beta = double( nrow_out * ncol_beta )
- # 'vcbeta' is square matrix with
- # width/height = (widthXL + widthXRs)
- # for now, will skip it
- # ...
- )$beta
- out <- matrix(out,ncol=ncol_beta)
- out
-}
\ No newline at end of file
Copied: pkg/MixABEL/tests/broken/rwth_caller.R (from rev 1958, pkg/MixABEL/R/rwth_caller.R)
===================================================================
--- pkg/MixABEL/tests/broken/rwth_caller.R (rev 0)
+++ pkg/MixABEL/tests/broken/rwth_caller.R 2015-05-10 07:44:20 UTC (rev 1967)
@@ -0,0 +1,60 @@
+"rwth_caller" <- function(
+ printall,
+# parameters and data to compute M
+ Phi, H2, Sigma2,
+# vector of length length_y
+ y,
+# matrix of height length_y and width widthXL
+ XL,
+# matrix of height length_y and width nXRs*widthXRs
+# it contains nXRs matrices XR
+ XRs,
+# width of single XR
+ WidthXR
+)
+{
+# figure out data dims and some sanity checks
+ Length_y = length(y);
+ NXRs = dim(XRs)[2]/WidthXR;
+ if ((dim(XRs)[2] %% WidthXR) != 0) stop("(dim(XRs)[2] %% widthXR) != 0");
+ if (dim(XRs)[1] != length(y)) stop("dim(XRs)[1] != length(y)");
+ if (is.null(XL)) {
+ XL <- matrix(1,ncol=1,nrow=length(y));
+ } else {
+ tmp <- matrix(1,ncol=1,nrow=length(y));
+ XL <- cbind(tmp,XL)
+ }
+ if (dim(XL)[1] != length(y)) stop("dim(XL)[1] != length(y)");
+ WidthXL = dim(XL)[2];
+
+ ncol_beta <- WidthXL + WidthXR;
+ nrow_out <- NXRs;
+ ncol_vcbeta <- ncol_beta*(ncol_beta+1)/2
+
+ out <- .C("rwth_example",
+ as.integer(printall),
+# dimensions of the data, explained later
+ as.integer(Length_y), as.integer(NXRs),
+ as.integer(WidthXR), as.integer(WidthXL),
+ # parameters and data to compute M
+ as.double(Phi), as.double(H2), as.double(Sigma2),
+ # vector of length length_y
+ as.double(y),
+ # matrix of height length_y and width widthXL
+ as.double(XL),
+ # matrix of height length_y and width nXRs*widthXRs
+ # it contains nXRs matrices XR
+ as.double(XRs),
+ # this is output
+ # 'beta's length is (widthXL + widthXRs)
+ # for time-being, just fill it with mean of
+ # columns of X
+ beta = double( nrow_out * ncol_beta )
+ # 'vcbeta' is square matrix with
+ # width/height = (widthXL + widthXRs)
+ # for now, will skip it
+ # ...
+ )$beta
+ out <- matrix(out,ncol=ncol_beta)
+ out
+}
\ No newline at end of file
Copied: pkg/MixABEL/tests/broken/test_rwth.R (from rev 1958, pkg/MixABEL/tests/test_rwth.R)
===================================================================
--- pkg/MixABEL/tests/broken/test_rwth.R (rev 0)
+++ pkg/MixABEL/tests/broken/test_rwth.R 2015-05-10 07:44:20 UTC (rev 1967)
@@ -0,0 +1,78 @@
+library(GenABEL)
+library(MixABEL)
+source("rwth_caller.R")
+
+data(ge03d2.clean)
+df <- ge03d2.clean[,autosomal(ge03d2.clean)]
+
+#load("saved_gkin_and_h2.RData")
+gkin <- ibs(df,w="freq")
+
+# toy example to see if the data transfer worked
+NIDS <- 7
+NSNPS <- 5
+phi <- matrix(c(1:(NIDS^2)),ncol=NIDS)
+phi
+XRs <- as.numeric(gtdata(df[1:NIDS,1:NSNPS]))
+XRs[is.na(XRs)] <- 0.5
+XRs
+y <- phdata(df)$height[1:NIDS]
+y[is.na(y)] <- mean(y,na.rm=T)
+y
+XL <- as.matrix(phdata(df)[1:NIDS,c("sex","age")])
+XL
+if (any(is.na(XL))) stop("opan'ki!")
+
+colMeans(XL)
+colMeans(XRs)
+a <- rwth_caller(1,phi,0.3,12.1,y,XL,XRs,1)
+a
+
+# repeat toy with different XR
+tmp <- matrix(NA,ncol=2*dim(XRs)[2],nrow=dim(XRs)[1])
+for (i in 1:dim(XRs)[2]) {
+ tmp[,i*2-1] <- 1*(XRs[,i]==1)
+ tmp[,i*2] <- 1*(XRs[,i]==2)
+}
+XRs <- tmp
+colMeans(XRs)
+a <- rwth_caller(1,phi,0.3,12.1,y,XL,XRs,2)
+a
+
+
+# realistic data
+# note that h2 is actually on lower end for this example
+# (where ~ 0.11, more frequently you see 0.25-0.75)
+# also phi is on 'sparse' end -- we may have even more sparse
+# data, but our typical phi is somewaht denser
+# (see 2,400 matrix Paolo has)
+# estimate model
+# h2ht <- polygenic(height~sex+age,data=df,kin=gkin)
+# save(gkin,h2ht,file="saved_gkin_and_h2.RData")
+# h2ht$h2an
+phi <- gkin[h2ht$measuredIDs,h2ht$measuredIDs]
+phi[upper.tri(phi)] <- t(phi)[upper.tri(phi)]
+phi[1:5,1:5]
+y <- h2ht$residualY[h2ht$measuredIDs]
+any(is.na(y))
+h2 <- h2ht$est
+h2
+Sigma2 <- var(y)
+Sigma2
+XL <- NULL
+XRs <- as.numeric(gtdata(df[h2ht$measuredIDs,1:10]))
+cmn <- colMeans(XRs,na.rm=T)
+for (i in 1:ncol(XRs)) {
+ XRs[is.na(XRs[,i]),i] <- cmn[i]
+}
+#XRs
+a <- rwth_caller(0,phi,h2,Sigma2,y,XL,XRs,1)
+cmn
+a
+
+old_res <- GWFGLS(y~SNP,data=phdata(df[h2ht$measuredIDs,]),
+ W=h2ht$InvSigma,genodata = XRs,verbosity=2,varcov=TRUE)
+# here beta_(Int.. and beta_SNP are elemnts of beta, and
+# se... make diagonal of 'e' and cov_.. is an off-diag elem
+# (so 'e' is 2x2 matrix but symmetric, so 3 numbers)
+old_res
Deleted: pkg/MixABEL/tests/test_rwth.R
===================================================================
--- pkg/MixABEL/tests/test_rwth.R 2015-05-08 13:00:25 UTC (rev 1966)
+++ pkg/MixABEL/tests/test_rwth.R 2015-05-10 07:44:20 UTC (rev 1967)
@@ -1,76 +0,0 @@
-library(GenABEL)
-library(MixABEL)
-data(ge03d2.clean)
-df <- ge03d2.clean[,autosomal(ge03d2.clean)]
-
-load("saved_gkin_and_h2.RData")
-#gkin <- ibs(df,w="freq")
-
-# toy example to see if the data transfer worked
-NIDS <- 7
-NSNPS <- 5
-phi <- matrix(c(1:(NIDS^2)),ncol=NIDS)
-phi
-XRs <- as.numeric(gtdata(df[1:NIDS,1:NSNPS]))
-XRs[is.na(XRs)] <- 0.5
-XRs
-y <- phdata(df)$height[1:NIDS]
-y[is.na(y)] <- mean(y,na.rm=T)
-y
-XL <- as.matrix(phdata(df)[1:NIDS,c("sex","age")])
-XL
-if (any(is.na(XL))) stop("opan'ki!")
-
-colMeans(XL)
-colMeans(XRs)
-a <- rwth_caller(1,phi,0.3,12.1,y,XL,XRs,1)
-a
-
-# repeat toy with different XR
-tmp <- matrix(NA,ncol=2*dim(XRs)[2],nrow=dim(XRs)[1])
-for (i in 1:dim(XRs)[2]) {
- tmp[,i*2-1] <- 1*(XRs[,i]==1)
- tmp[,i*2] <- 1*(XRs[,i]==2)
-}
-XRs <- tmp
-colMeans(XRs)
-a <- rwth_caller(1,phi,0.3,12.1,y,XL,XRs,2)
-a
-
-
-# realistic data
-# note that h2 is actually on lower end for this example
-# (where ~ 0.11, more frequently you see 0.25-0.75)
-# also phi is on 'sparse' end -- we may have even more sparse
-# data, but our typical phi is somewaht denser
-# (see 2,400 matrix Paolo has)
-# estimate model
-# h2ht <- polygenic(height~sex+age,data=df,kin=gkin)
-# save(gkin,h2ht,file="saved_gkin_and_h2.RData")
-# h2ht$h2an
-phi <- gkin[h2ht$measuredIDs,h2ht$measuredIDs]
-phi[upper.tri(phi)] <- t(phi)[upper.tri(phi)]
-phi[1:5,1:5]
-y <- h2ht$residualY[h2ht$measuredIDs]
-any(is.na(y))
-h2 <- h2ht$est
-h2
-Sigma2 <- var(y)
-Sigma2
-XL <- NULL
-XRs <- as.numeric(gtdata(df[h2ht$measuredIDs,1:10]))
-cmn <- colMeans(XRs,na.rm=T)
-for (i in 1:ncol(XRs)) {
- XRs[is.na(XRs[,i]),i] <- cmn[i]
-}
-#XRs
-a <- rwth_caller(0,phi,h2,Sigma2,y,XL,XRs,1)
-cmn
-a
-
-old_res <- GWFGLS(y~SNP,data=phdata(df[h2ht$measuredIDs,]),
- W=h2ht$InvSigma,genodata = XRs,verbosity=2,varcov=TRUE)
-# here beta_(Int.. and beta_SNP are elemnts of beta, and
-# se... make diagonal of 'e' and cov_.. is an off-diag elem
-# (so 'e' is 2x2 matrix but symmetric, so 3 numbers)
-old_res
More information about the Genabel-commits
mailing list