### R code from vignette source 'vignettes/globaltest/inst/doc/GlobalTest.Rnw' ################################################### ### code chunk number 1: make data ################################################### set.seed(1) Y <- rnorm(20) X <- matrix(rnorm(200), 20, 10) X[,1:3] <- X[,1:3] + Y colnames(X) <- LETTERS[1:10] ################################################### ### code chunk number 2: GlobalTest.Rnw:98-99 ################################################### library(globaltest) ################################################### ### code chunk number 3: options ################################################### gt.options(trace=FALSE, max.print=45) ################################################### ### code chunk number 4: full formula simple ################################################### gt(Y~1, Y~A+B+C, data = X) ################################################### ### code chunk number 5: full formula ################################################### gt(Y~D, Y~A+B+C+D, data = X) ################################################### ### code chunk number 6: summary ################################################### summary(gt(Y~A, Y~A+B+C, data = X)) ################################################### ### code chunk number 7: GlobalTest.Rnw:140-145 ################################################### res <- gt(Y~A, Y~A+B+C, data = X) p.value(res) z.score(res) result(res) size(res) ################################################### ### code chunk number 8: GlobalTest.Rnw:156-157 (eval = FALSE) ################################################### ## gt(Y~A, ~B+C, data = X) ################################################### ### code chunk number 9: GlobalTest.Rnw:162-163 ################################################### gt(Y~A+B+C, data = X) ################################################### ### code chunk number 10: GlobalTest.Rnw:167-168 ################################################### gt(Y~A, ~., data = X) ################################################### ### code chunk number 11: GlobalTest.Rnw:173-174 ################################################### gt(Y, X) ################################################### ### code chunk number 12: null design ################################################### designA <- cbind(1, X[,"A"]) gt(Y, X, designA) ################################################### ### code chunk number 13: poisson ################################################### P <- rpois(20, lambda=2) gt(P~A, ~., data=X, model = "Poisson") gt(P~A, ~., data=X, model = "linear") ################################################### ### code chunk number 14: permutations ################################################### gt(Y,X) gt(Y,X, permutations=1e4) ################################################### ### code chunk number 15: GlobalTest.Rnw:221-222 ################################################### hist(gt(Y,X, permutations=1e4)) ################################################### ### code chunk number 16: GlobalTest.Rnw:231-234 ################################################### A <- X[,"A"] gt(Y,X,A) gt(Y,X,~A) ################################################### ### code chunk number 17: no intercept ################################################### gt(Y~0+A, ~ B+C, data = X) ################################################### ### code chunk number 18: alternative intercept ################################################### IC <- rep(1, 20) gt(Y~0+A, ~ IC+B+C, data = X) ################################################### ### code chunk number 19: factors ################################################### YY <- rnorm(6) FF <- factor(rep(letters[1:2], 3)) GG <- factor(rep(letters[3:5], 2)) model.matrix(gt(YY ~ FF + GG, x = TRUE))$alternative ################################################### ### code chunk number 20: ordered factors ################################################### GG <- ordered(GG) model.matrix(gt(YY ~ GG, x = TRUE))$alternative ################################################### ### code chunk number 21: weights ################################################### res <- gt(Y, X) weights(res) ################################################### ### code chunk number 22: GlobalTest.Rnw:281-283 ################################################### res <- gt(Y,X, standardize=TRUE) weights(res) ################################################### ### code chunk number 23: GlobalTest.Rnw:287-289 (eval = FALSE) ################################################### ## gt(Y, X[,c("A","A","B")], weights=c(.5,.5,1)) ## gt(Y, X[,c("A","B")]) ################################################### ### code chunk number 24: GlobalTest.Rnw:296-297 ################################################### gt(Y, X, directional = TRUE) ################################################### ### code chunk number 25: GlobalTest.Rnw:308-309 ################################################### gt(Y~A+B+C,data=X, test.value=c(.2,.2,.2)) ################################################### ### code chunk number 26: GlobalTest.Rnw:314-316 ################################################### os <- X[,1:3]%*%c(.2,.2,.2) gt(Y~offset(os), ~A+B+C, data=X) ################################################### ### code chunk number 27: covariates preparation ################################################### gt(Y~A+B, data=X) gt(Y~A, data=X) gt(Y~B, data=X) ################################################### ### code chunk number 28: covariates ################################################### covariates(gt(Y,X)) ################################################### ### code chunk number 29: covariatesW ################################################### covariates(gt(Y,X), what="w") ################################################### ### code chunk number 30: covariates output ################################################### res <- covariates(gt(Y,X)) res[1:10] ################################################### ### code chunk number 31: covariates leafNodes ################################################### leafNodes(res, alpha=0.10) ################################################### ### code chunk number 32: extract ################################################### extract(res) ################################################### ### code chunk number 33: covariates_zoom ################################################### covariates(gt(Y,X), zoom=TRUE) ################################################### ### code chunk number 34: GlobalTest.Rnw:405-406 ################################################### subjects(gt(Y,X)) ################################################### ### code chunk number 35: GlobalTest.Rnw:413-414 ################################################### subjects(gt(Y,X), what="s", mirror=FALSE) ################################################### ### code chunk number 36: subset ################################################### set <- LETTERS[1:3] gt(Y,X, subsets = set) ################################################### ### code chunk number 37: many sets ################################################### sets <- list(one=LETTERS[1:3], two=LETTERS[4:6]) gt(Y,X, subsets = sets) ################################################### ### code chunk number 38: subsets method ################################################### res <- gt(Y,X, subsets = sets) subsets(res) ################################################### ### code chunk number 39: many weights ################################################### wts <- list(up = 1:10, down = 10:1) gt(Y,X,weights=wts) ################################################### ### code chunk number 40: single weight many subsets ################################################### gt(Y,X, subsets=sets, weights=1:10) ################################################### ### code chunk number 41: subsets and weights ################################################### gt(Y,X, subsets=sets, weights=wts) gt(Y,X, subsets=sets, weights=list(1:3,7:5)) ################################################### ### code chunk number 42: alias ################################################### res <- gt(Y,X, weights=wts, alias = c("one", "two")) alias(res) alias(res) <- c("ONE", "TWO") ################################################### ### code chunk number 43: sort ################################################### res[1] sort(res) ################################################### ### code chunk number 44: p.adjust ################################################### p.adjust(res) p.adjust(res, "BH") p.adjust(res, "BY") ################################################### ### code chunk number 45: focuslevel set structure ################################################### level1 <- as.list(LETTERS[1:10]) names(level1) <- letters[1:10] level2 <- list(abc = LETTERS[1:3], cde = LETTERS[3:5], fgh = LETTERS[6:8], hij = LETTERS[8:10]) level3 <- list(all = LETTERS[1:10]) dag <- c(level1, level2, level3) ################################################### ### code chunk number 46: focus level choice ################################################### fl <- names(level2) fl <- findFocus(dag, maxsize=8) ################################################### ### code chunk number 47: focus level ################################################### res <- gt(Y,X) res <- focusLevel(res, sets = dag, focus=fl) sort(res) ################################################### ### code chunk number 48: focus level leaf nodes ################################################### leafNodes(res) ################################################### ### code chunk number 49: draw ################################################### draw(res, names=TRUE) ################################################### ### code chunk number 50: draw legend (eval = FALSE) ################################################### ## legend <- draw(res) ################################################### ### code chunk number 51: inheritance set structure ################################################### level1 <- as.list(LETTERS[1:10]) names(level1) <- letters[1:10] level2 <- list(ab = LETTERS[1:2], cde = LETTERS[3:5], fg = LETTERS[6:7], hij = LETTERS[8:10]) level3 <- list(all = LETTERS[1:10]) tree <- c(level1, level2, level3) ################################################### ### code chunk number 52: inheritance ################################################### res <- gt(Y,X) resI <- inheritance(res, tree) resI ################################################### ### code chunk number 53: inheritance ################################################### hc <- hclust(dist(t(X))) resHC <- inheritance(res, hc) resHC ################################################### ### code chunk number 54: inheritance leaf nodes ################################################### leafNodes(resI) leafNodes(resHC) ################################################### ### code chunk number 55: inheritance_draw ################################################### draw(resHC, names=TRUE) ################################################### ### code chunk number 56: inheritance covariates (eval = FALSE) ################################################### ## covariates(res) ################################################### ### code chunk number 57: transpose option (eval = FALSE) ################################################### ## gt.options(transpose=TRUE) ################################################### ### code chunk number 58: load Golub ################################################### library(golubEsets) data(Golub_Train) ################################################### ### code chunk number 59: vsn ################################################### library(vsn) exprs(Golub_Train) <- exprs(vsn2(Golub_Train)) ################################################### ### code chunk number 60: overall ################################################### gt(ALL.AML, Golub_Train) ################################################### ### code chunk number 61: Source ################################################### gt(ALL.AML ~ Source, Golub_Train) ################################################### ### code chunk number 62: transpose option (eval = FALSE) ################################################### ## gt.options(transpose=TRUE) ################################################### ### code chunk number 63: trim option (eval = FALSE) ################################################### ## gt.options(trim=TRUE) ################################################### ### code chunk number 64: gtKEGG ################################################### gtKEGG(ALL.AML, Golub_Train, id = "04110") ################################################### ### code chunk number 65: KEGG organism package ################################################### eg <- as.list(hu6800ENTREZID) gtKEGG(ALL.AML, Golub_Train, id="04110", probe2entrez = eg, annotation="org.Hs.eg.db") ################################################### ### code chunk number 66: gtKEGG multtest ################################################### gtKEGG(ALL.AML, Golub_Train, id=c("04110","04210"), multtest="BH") ################################################### ### code chunk number 67: testKEGG (eval = FALSE) ################################################### ## gtKEGG(ALL.AML, Golub_Train) ################################################### ### code chunk number 68: gtGO ################################################### gtGO(ALL.AML, Golub_Train, id="GO:0007049") ################################################### ### code chunk number 69: GO organism package ################################################### eg <- as.list(hu6800ENTREZID) gtGO(ALL.AML, Golub_Train, id="GO:0007049", probe2entrez = eg, annotation="org.Hs.eg") ################################################### ### code chunk number 70: testBP (eval = FALSE) ################################################### ## gtGO(ALL.AML, Golub_Train, ontology="BP", minsize = 10, maxsize = 500) ################################################### ### code chunk number 71: gtGO multtest ################################################### gtGO(ALL.AML, Golub_Train, id=c("GO:0007049","GO:0006915"), multtest="BH") ################################################### ### code chunk number 72: focusGO ################################################### descendants <- get("GO:0007049", GOBPOFFSPRING) res <- gtGO(ALL.AML, Golub_Train, id = c("GO:0007049", descendants), multtest = "focus") leafNodes(res) ################################################### ### code chunk number 73: significantGO1 (eval = FALSE) ################################################### ## draw(res, interactive=TRUE) ## legend <- draw(res) ################################################### ### code chunk number 74: significantGO2 ################################################### draw(res) ################################################### ### code chunk number 75: getBroadSets (eval = FALSE) ################################################### ## broad <- getBroadSets("your/path/to/msigdb_v.2.5.xml") ################################################### ### code chunk number 76: gtBroad (eval = FALSE) ################################################### ## gtBroad(ALL.AML, Golub_Train, id = "chr5q33", collection=broad) ################################################### ### code chunk number 77: Broad organism package (eval = FALSE) ################################################### ## eg <- as.list(hu6800ENTREZID) ## gtBroad(ALL.AML, Golub_Train, id = "chr5q33", collection=broad, probe2entrez = eg, annotation="org.Hs.eg.db") ################################################### ### code chunk number 78: gtBroad multtest (eval = FALSE) ################################################### ## gtBroad(ALL.AML, Golub_Train, id=c("chr5q33","chr5q34"), multtest="BH", collection=broad) ################################################### ### code chunk number 79: Broad c1 (eval = FALSE) ################################################### ## gtBroad(ALL.AML, Golub_Train, category="c1", collection=broad) ################################################### ### code chunk number 80: gtConcept (eval = FALSE) ################################################### ## gtConcept(ALL.AML, Golub_Train, conceptmatrix="Body System.txt") ################################################### ### code chunk number 81: gtConcept organism package (eval = FALSE) ################################################### ## eg <- as.list(hu6800ENTREZID) ## gtConcept(ALL.AML, Golub_Train, conceptmatrix="Body System.txt", probe2entrez = eg, annotation="org.Hs.eg.db") ################################################### ### code chunk number 82: gtConcept multtest (eval = FALSE) ################################################### ## gtConcept(ALL.AML, Golub_Train, conceptmatrix="Body System.txt", multtest="BH") ################################################### ### code chunk number 83: KEGG_covariates ################################################### res <- gtKEGG(ALL.AML, Golub_Train, id = "04110") features(res, alias=hu6800SYMBOL) ################################################### ### code chunk number 84: LeafNodes geneplot ################################################### ft <- features(res, alias=hu6800SYMBOL) leafNodes(ft) ################################################### ### code chunk number 85: subsets leafNodes (eval = FALSE) ################################################### ## subsets(leafNodes(ft)) ################################################### ### code chunk number 86: KEGG_covariates_zoom ################################################### res <- gtKEGG(ALL.AML, Golub_Train, id = "04110") features(res, alias=hu6800SYMBOL, zoom=TRUE) ################################################### ### code chunk number 87: extract geneset ################################################### ft <- features(res, alias=hu6800SYMBOL, zoom=TRUE) extract(ft) ################################################### ### code chunk number 88: pdf covariates (eval = FALSE) ################################################### ## res_all <- gtKEGG(ALL.AML, Golub_Train) ## features(res_all[1:5], pdf="KEGGcov.pdf", alias=hu6800SYMBOL) ################################################### ### code chunk number 89: KEGG_subjects ################################################### res <- gt(ALL.AML, Golub_Train) subjects(res) ################################################### ### code chunk number 90: pdf covariates (eval = FALSE) ################################################### ## res_all <- gtKEGG(ALL.AML, Golub_Train) ## subjects(res_all[1:25], pdf="KEGGsubj.pdf") ################################################### ### code chunk number 91: lungExpression ################################################### library(lungExpression) data(michigan) gt(Surv(TIME..months., death==1), michigan) ################################################### ### code chunk number 92: comparative ################################################### res <- gtKEGG(ALL.AML, Golub_Train, id = "04110") comparative(res) ################################################### ### code chunk number 93: faults ################################################### require("boot") data(cloth) Z <- matrix(diag(nrow(cloth)), ncol = nrow(cloth), dimnames = list(NULL, 1:nrow(cloth))) gt(y ~ log(x), alternative = Z, data = cloth, model = "poisson") ################################################### ### code chunk number 94: litters ################################################### library("survival") data(rats) nlitters<-length(unique(rats$litter)) Z<-matrix(NA,dim(rats)[1],nlitters, dimnames=list(NULL,1:nlitters)) for (i in 1:nlitters) Z[,i]<-(rats[,1]==i)*1 gt(Surv(time,status)~rx,alternative=Z,data=rats,model="cox") ################################################### ### code chunk number 95: anscombe ################################################### data(anscombe) set.seed(0) X<-anscombe$x2 Y<-anscombe$y2 + rnorm(length(X),0,3) gtPS(Y~X) ################################################### ### code chunk number 96: reparamZ ################################################### Z<-bbase(X,bdeg=3,nint=10) P<-reparamZ(Z,pord=2) gt(Y~X,alternative=P) ################################################### ### code chunk number 97: Uspan ################################################### U<-reparamZ(Z,pord=2, returnU=TRUE)$U lm(Y~X+U)$coefficients ################################################### ### code chunk number 98: arg pord ################################################### gtPS(Y~X, pord=list(Z=0,P=2)) ################################################### ### code chunk number 99: arg robust ################################################### rob<-gtPS(Y~X, pord=list(Z=0,P=2), robust=TRUE) rob@result ################################################### ### code chunk number 100: cbind weights ################################################### comb<-gt(Y~X, alternative=cbind(Z,P)) comb@result ################################################### ### code chunk number 101: robust weights ################################################### colrange<-list(Z=1:ncol(Z), P=(ncol(Z)+1):(ncol(Z)+ncol(P))) sapply(list(combined=comb,robust=rob), function(x){sapply(colrange, function(y){sum(weights(x)[y])/sum(weights(x))})}) ################################################### ### code chunk number 102: kyphosis ################################################### require("rpart") data("kyphosis") fit0<-glm(Kyphosis~., data = kyphosis, family="binomial") res<-gtPS(fit0) res@result ################################################### ### code chunk number 103: sterms ################################################### sterms(res) ################################################### ### code chunk number 104: penalized fit ################################################### require("penalized") Z<-gtPS(fit0, returnZ=TRUE)$Z fit1<-penalized(Kyphosis, penalized=~ Z, unpenalized=~Age+Number+Start, data = kyphosis, model="logistic", lambda2 = 0.086, trace=FALSE) ################################################### ### code chunk number 105: goodness_kyphosis ################################################### bd=3 ni=10 po=2 covnames<-names(kyphosis)[-1] d<-bd+ni-po gammas<-fit1@penalized betas<-fit1@unpenalized l<-length(covnames) cd<-c(0,cumsum(rep(d,l))) op <- par(mfrow = c(2, 2)) for (i in 1:3){ x<-kyphosis[,covnames[i]] sx<-sort(x,index.return=T) ind<-vector("list", l) ind[[i]]<-(cd[i]+1):(cd[i+1]) plot(sx$x, (Z[sx$ix,ind[[i]]]%*%gammas[ind[[i]]] ), type="b",ylim=c(-1.7,0.7), xlab=covnames[i], ylab=paste("s(",covnames[i],")", sep="")) lines(abline(h=0,lty="dotted")) rug(x) } par(op) ################################################### ### code chunk number 106: arg covs ################################################### gtPS(fit0, covs=c("Age","Start")) ################################################### ### code chunk number 107: arg nint ################################################### gtPS(fit0,covs=c("Age","Number","Start"), nint=list(a=5, b=c(5,1,1)), pord=0) ################################################### ### code chunk number 108: cbind GAM ################################################### covs=c("Age","Number","Start") bd=c(3,3,3);ni=c(10,10,10);po=c(2,2,2);cs<-c(0,cumsum(bd+ni-po)) X0<-model.matrix(fit0)[,] combZ<-do.call(cbind,lapply(1:length(covs),function(x){reparamZ(bbase(kyphosis[,covs[x]], nint=ni[x], bdeg=bd[x]), pord=po[x])})) comb<-gt(Kyphosis~., alternative=combZ, data = kyphosis, model="logistic") comb@result ################################################### ### code chunk number 109: cbind weights GAM ################################################### range<-lapply(1:length(covs),function(x){(cs[x]+1):(cs[x+1])}) names(range)<-covs sapply(range,function(x){sum(weights(comb)[x])/sum(weights(comb))}) ################################################### ### code chunk number 110: reweighZ ################################################### rwgtZ<-do.call(cbind,lapply(1:length(covs),function(x){reweighZ(reparamZ(bbase(kyphosis[,covs[x]], nint=ni[x], bdeg=bd[x]), pord=po[x]),fit0)})) rwgt<-gt(Kyphosis~., alternative=rwgtZ, data = kyphosis, model="logistic") sapply(range,function(x){sum(weights(rwgt)[x])/sum(weights(rwgt))}) ################################################### ### code chunk number 111: lake ################################################### library(gss) data(LakeAcidity) fit0<-lm(ph~log10(cal)+lat+lon, data=LakeAcidity) res<-gtKS(fit0) res@result sterms(res) ################################################### ### code chunk number 112: trace ################################################### gtKS(fit0, quant=seq(.01,.99,.02), data=LakeAcidity, termlabels=TRUE, robust=T) ################################################### ### code chunk number 113: goodness_trace ################################################### p<-sapply(seq(.01,.99,.02),function(x){gtKS(fit0, termlabels=T,quant=x, data=LakeAcidity)@result[,1]}) plot(seq(.01,.99,.02),p, type="s", xlab="quant", ylab="p-value", ylim=c(0,.5), xlim=c(0,1)) abline(h=gtKS(ph~log10(cal)+lat+lon, quant=seq(.01,.99,.02), data=LakeAcidity, robust=T, termlabels=T)@result[,1], lty="dotted") ################################################### ### code chunk number 114: btensor ################################################### fit0<-lm(ph~lat*lon, data=LakeAcidity) res<-gtPS(fit0, covs=c("lat","lon"), interact=TRUE, data=LakeAcidity) res@result sterms(res) ################################################### ### code chunk number 115: goodness_lake ################################################### Z<-gtPS(fit0, returnZ=TRUE, interact=TRUE)$Z fit1<-penalized(ph, penalized=Z, unpenalized=~lat*lon, data = LakeAcidity, lambda2 = 1, trace=FALSE) lon<-seq(79,85,length=50) lat<-seq(33,38,length=50) new<-data.frame(expand.grid(lon=lon, lat=lat)) Znew<-btensor(new, bdeg=c(3,3), nint=c(10,10), pord=c(2,2)) fitted<-matrix(predict(fit1,unpenalized=~lat*lon, penalized=Znew, data=new)[,1],50,50) persp(lon,lat,fitted,theta=-25) ################################################### ### code chunk number 116: anova type ################################################### Z1<-reweighZ(reparamZ(bbase(LakeAcidity$lat, bdeg=3, nint=10), pord=2), fit0) Z2<-reweighZ(reparamZ(bbase(LakeAcidity$lon, bdeg=3, nint=10), pord=2), fit0) Z12<-reweighZ(btensor(cbind(LakeAcidity$lat, LakeAcidity$lon),bdeg=c(3,3),nint=c(10,10),pord=c(2,2)), fit0) gt(ph~lat*lon, alternative=cbind(Z1,Z2,Z12), data=LakeAcidity) ################################################### ### code chunk number 117: varying ################################################### data(nox) sE<-bbase(nox$equi, bdeg=3, nint=10) sEbyC<-model.matrix(~0+sE:factor(comp), data=nox)[,] gt(nox~equi*factor(comp), alternative=cbind(sE,sEbyC), data=nox) ################################################### ### code chunk number 118: boston ################################################### library(MASS) data(Boston) res<-gtLI(medv~., data=Boston) res@result round(weights(res)/sum(weights(res)),4) ################################################### ### code chunk number 119: GlobalTest.Rnw:1242-1244 ################################################### gtLI(medv~., data=Boston, standardize=T) gtLI(medv~., data=scale(Boston))