## ----include=FALSE------------------------------------------------------------ knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----eval=FALSE, include=FALSE------------------------------------------------ # install.packages("itdr") # library(itdr) ## ----eval=FALSE, include=TRUE------------------------------------------------- # install.packages("~/itdr.zip") # library(itdr) ## ----eval=FALSE, include=TRUE------------------------------------------------- # library(devtools) # install_github("TharinduPDeAlwis/itdr") # library(itdr) ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data(automobile) automobile.na <- na.omit(automobile) # prepare response and predictor variables auto_y <- log(automobile.na[, 26]) auto_xx <- automobile.na[, c(10, 11, 12, 13, 14, 17, 19, 20, 21, 22, 23, 24, 25)] auto_x <- scale(auto_xx) # Standardize the predictors # call to the d.boots() function with required #arguments d_est <- d.boots(auto_y, auto_x, var_plot = TRUE, space = "pdf", xdensity = "normal", method = "FM" ) auto_d <- d_est$d.hat auto_d ## ----eval=TRUE, include=TRUE-------------------------------------------------- auto_d <- 2 # The estimated value from Section 2.1 candidate_list<-seq(0.05, 1, by = 0.01) auto_sw2 <- hyperPara(auto_y, auto_x, auto_d, range = candidate_list,xdensity = "normal", B = 500, space = "pdf", method = "FM", hyper = "wx") auto_sw2$wx.hat # we get the estimator for sw2 as 0.09 ## ----eval=TRUE, include=TRUE-------------------------------------------------- auto_d <- 2 # Estimated value from Section 2.1 auto_st2 <- hyperPara(auto_y, auto_x, auto_d, wx = 0.1, range = seq(0.1, 1, by = 0.1), xdensity = "normal", method = "FM", hyper = "wy") auto_st2$wy.hat # we get the estimator for st2=0.9 ## ----eval=TRUE, include=TRUE-------------------------------------------------- h_hat <- hyperPara(auto_y, auto_x, auto_d, wx = 5, wy = 0.1, range = seq(0.1, 2, by = .1), space = "pdf", method = "FM", hyper = "wh") # Bandwidth estimator for Gaussian kernel density estimation for central subspace h_hat$h.hat # we have the estimator as h_hat=0.1 ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data(automobile) head(automobile) df <- cbind(automobile[, c(26, 10, 11, 12, 13, 14, 17, 19, 20, 21, 22, 23, 24, 25)]) dff <- as.matrix(df) automobi <- dff[complete.cases(dff), ] d <- 2 # Estimated value from Section 2.1 wx <- .14 # Estimated value from Section 2.2.1 wy <- .9 # Estimated value from Section 2.2.2 wh <- 1.5 # Estimated value from Section 2.2.3 p <- 13 # Estimated value from Section 2.3 y <- automobi[, 1] x <- automobi[, c(2:14)] xt <- scale(x) # Distribution of the predictors is a normal distribution fit.F_CMS <- itdr(y, xt, d, wx, wy, wh, space = "pdf", xdensity = "normal", method = "FM") round(fit.F_CMS$eta_hat, 2) # Distribution of the predictors is a unknown (using kernel method) fit.F_CMS <- itdr(y, xt, d, wx, wy, wh, space = "pdf", xdensity = "kernel", method = "FM") round(fit.F_CMS$eta_hat, 2) ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data("recumbent") recumbent.na <- na.omit(recumbent) y <- recumbent.na$outcome X1 <- log(recumbent.na$ast) X2 <- log(recumbent.na$ck) X3 <- log(recumbent.na$urea) p <- 3 x <- matrix(c(X1, X2, X3), ncol = p) d <- 2 fit.iht_CMS <- itdr(y, x, 2, method = "iht") fit.iht_CMS$eta_hat ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data(pdb) colnames(pdb) <- NULL p <- 15 # select predictor vecotr (y) and response variables (X) according to Weng and Weng and Yin, (2018). df <- pdb[, c(79, 73, 77, 103, 112, 115, 124, 130, 132, 145, 149, 151, 153, 155, 167, 169)] dff <- as.matrix(df) # remove the NA rows planingdb <- dff[complete.cases(dff), ] y <- planingdb[, 1] # n-dimensionl response vector x <- planingdb[, c(2:(p + 1))] # raw desing matrix x <- x + 0.5 # desing matrix after tranformations xt <- cbind( x[, 1]^(.33), x[, 2]^(.33), x[, 3]^(.57), x[, 4]^(.33), x[, 5]^(.4), x[, 6]^(.5), x[, 7]^(.33), x[, 8]^(.16), x[, 9]^(.27), x[, 10]^(.5), x[, 11]^(.5), x[, 12]^(.33), x[, 13]^(.06), x[, 14]^(.15), x[, 15]^(.1) ) m <- 1 W <- sapply(50, rnorm) # run the hypothsis tests d.test(y, x, m) ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data(pdb) colnames(pdb) <- NULL p <- 15 # select predictor vecotr (y) and response variables (X) according to Weng and Weng and Yin, (2018). df <- pdb[, c(79, 73, 77, 103, 112, 115, 124, 130, 132, 145, 149, 151, 153, 155, 167, 169)] dff <- as.matrix(df) # remove the NA rows planingdb <- dff[complete.cases(dff), ] y <- planingdb[, 1] # n-dimensionl response vector x <- planingdb[, c(2:(p + 1))] # raw desing matrix x <- x + 0.5 # desing matrix after tranformations give in Weng and Yin, (2018). xt <- cbind( x[, 1]^(.33), x[, 2]^(.33), x[, 3]^(.57), x[, 4]^(.33), x[, 5]^(.4), x[, 6]^(.5), x[, 7]^(.33), x[, 8]^(.16), x[, 9]^(.27), x[, 10]^(.5), x[, 11]^(.5), x[, 12]^(.33), x[, 13]^(.06), x[, 14]^(.15), x[, 15]^(.1) ) d <- 1 # estimated dimension of the CS from Section 4.1 invFM.fit <- itdr(y, x, d, m = 50, method = "invFM", x.scale = FALSE) # estimated basis betahat <- invFM.fit$eta_hat plot(y ~ xt %*% betahat, xlab = "First reduced predictor", ylab = "Health insurance coverage" ) ## ----eval=TRUE, include=TRUE-------------------------------------------------- library(itdr) data(prostate) X <- as.matrix(prostate[, 1:8]) Y <- matrix(prostate[, 9], ncol = 1) fit.ftire <- mitdr(X, Y, d = 2, m = 10, method = "FT-IRE") betahat <- fit.ftire$Beta_hat betahat newx <- X %*% betahat plot(Y ~ newx[, 1], xlab = "First reduced predictor", ylab = paste0(expression(log), "(antigen)", sep = "") ) plot(Y ~ newx[, 2], xlab = "Second reduced predictor", ylab = paste0(expression(log), "(antigen)", sep = "") ) ## ----eval=TRUE, include=TRUE-------------------------------------------------- data(raman) Y <- as.matrix(raman[, c(1100)]) ## percentage of total fat content X <- as.matrix(raman[c(2:501)]) ## first 500 wavelength variables out <- mitdr(X, Y, d = 1, m = 30, method = "admmft", lambda = 0.5, sparse.cov = TRUE, x.scale = TRUE) estbeta <- out$Beta_hat estbeta plot(Y ~ X %*% estbeta, xlab = "First reduced predictor", ylab = "Percentage of total fat" )