11
2- R version 4.0.3 (2020-10-10 ) -- "Bunny-Wunnies Freak Out "
3- Copyright (C) 2020 The R Foundation for Statistical Computing
2+ R Under development (unstable) (2021-04-05 r80145 ) -- "Unsuffered Consequences "
3+ Copyright (C) 2021 The R Foundation for Statistical Computing
44Platform: x86_64-pc-linux-gnu (64-bit)
55
66R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -1345,7 +1345,7 @@ In bols(z2, by = z3) : Dropped unobserved factor levels
13451345+ unzip("bkernel.zip")
13461346+ txt <- readLines("Kernel_Boosting_example_code.R")
13471347+ writeLines(txt[-c(1:10, 149:length(txt))], con = "run.R")
1348- + source("run.R", echo = TRUE )
1348+ + source("run.R", echo = FALSE )
13491349+ setwd(wd)
13501350+ }
13511351Loading required package: kangar00
@@ -1356,17 +1356,6 @@ The following object is masked from 'asNamespace("mboost")':
13561356
13571357 make_psd
13581358
1359-
1360- > set.seed(1104)
1361-
1362- > CORES <- 1
1363-
1364- > download <- FALSE
1365-
1366- > if (!require("KEGGgraph")) {
1367- + source("https://bioconductor.org/biocLite.R")
1368- + biocLite("KEGGgraph", suppressUpdates = TRUE)
1369- + }
13701359Loading required package: KEGGgraph
13711360
13721361Attaching package: 'KEGGgraph'
@@ -1379,113 +1368,21 @@ The following object is masked from 'package:base':
13791368
13801369 plot
13811370
1382-
1383- > if (!require("biomaRt")) {
1384- + source("https://bioconductor.org/biocLite.R")
1385- + biocLite("biomaRt", suppressUpdates = TRUE)
1386- + }
13871371Loading required package: biomaRt
1388-
1389- > if (!require("mboost")) {
1390- + install.packages("mboost")
1391- + }
1392-
1393- > if (!require("kangar00")) {
1394- + install.packages("kangar00")
1395- + }
1396-
1397- > library("mboost")
1398-
1399- > library("kangar00")
1400-
1401- > if (download == TRUE) {
1402- + list.of.pathways <- list()
1403- + for (pw in c("hsa04662", "hsa04917", "hsa05145", "hsa05210",
1404- + "hsa04024")) .... [TRUNCATED]
1405-
1406- > anno <- get_anno(snp.info, pathway.info)
1407-
1408- > geno <- get(load("genotype.data.Rda"))
1409-
1410- > pheno <- get(load("phenotype.data.Rda"))
1411-
1412- > gwd <- GWASdata(pheno = pheno, anno = anno, geno = geno,
1413- + desc = "study_data")
1414-
1415- > data <- list(case_control_status = pheno$case_control_status,
1416- + gwd = gwd)
1417-
1418- > kernels <- paste0("bkernel(gwd, kernel = \"net\", pathway = list.of.pathways[[",
1419- + seq(1:length(list.of.pathways)), "]], df = 4)")
1420-
1421- > kernels <- paste(kernels, collapse = " + ")
1422-
1423- > (fm <- as.formula(paste("case_control_status ~ ",
1424- + kernels)))
1425- case_control_status ~ bkernel(gwd, kernel = "net", pathway = list.of.pathways[[1]],
1426- df = 4) + bkernel(gwd, kernel = "net", pathway = list.of.pathways[[2]],
1427- df = 4) + bkernel(gwd, kernel = "net", pathway = list.of.pathways[[3]],
1428- df = 4) + bkernel(gwd, kernel = "net", pathway = list.of.pathways[[4]],
1429- df = 4) + bkernel(gwd, kernel = "net", pathway = list.of.pathways[[5]],
1430- df = 4)
1431-
1432- > data$case_control_status <- factor(data$case_control_status)
1433-
1434- > logit <- glm(case_control_status ~ age + sex, data = gwd@pheno,
1435- + family = binomial())
1436-
1437- > summary(logit)
1438-
1439- Call:
1440- glm(formula = case_control_status ~ age + sex, family = binomial(),
1441- data = gwd@pheno)
1442-
1443- Deviance Residuals:
1444- Min 1Q Median 3Q Max
1445- -1.40863 -1.16853 0.00252 1.15372 1.42560
1446-
1447- Coefficients:
1448- Estimate Std. Error z value Pr(>|z|)
1449- (Intercept) 0.800774 0.464828 1.723 0.0849 .
1450- age -0.012361 0.008824 -1.401 0.1612
1451- sex -0.452833 0.286332 -1.581 0.1138
1452- ---
1453- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1454-
1455- (Dispersion parameter for binomial family taken to be 1)
1456-
1457- Null deviance: 277.26 on 199 degrees of freedom
1458- Residual deviance: 272.85 on 197 degrees of freedom
1459- AIC: 278.85
1460-
1461- Number of Fisher Scoring iterations: 4
1462-
1463-
1464- > mod <- mboost(fm, data = data, control = boost_control(mstop = 0,
1465- + trace = TRUE), offset = predict(logit), family = Binomial())
1466-
1467- > cat("Cross-validation started.\n")
14681372Cross-validation started.
1469-
1470- > folds_subsampl <- cv(model.weights(mod), type = "subsampling",
1471- + B = 5)
1472-
1473- > cvr <- cvrisk(mod, folds = folds_subsampl, grid = 0:20,
1474- + mc.cores = CORES)
14751373[ 1] ..................
1476- Final risk: 68.80205
1374+ Final risk: 74.92605
14771375[ 1] ..................
1478- Final risk: 72.01955
1376+ Final risk: 66.94996
14791377[ 1] ..................
1480- Final risk: 67.02656
1378+ Final risk: 69.67069
14811379[ 1] ..................
1482- Final risk: 69.86766
1380+ Final risk: 66.76257
14831381[ 1] ..................
1484- Final risk: 71.82216
1485-
1486- > plot(cvr)
1487- There were 22 warnings (use warnings() to see them)
1382+ Final risk: 70.07644
1383+ Warning message:
1384+ no DISPLAY variable so Tk is not available
14881385>
14891386> proc.time()
14901387 user system elapsed
1491- 91.431 1.478 94.153
1388+ 96.621 2.589 101.235
0 commit comments