diff --git a/NAMESPACE b/NAMESPACE
index ff0401e..6349f00 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,6 +14,7 @@ export(make_group)
 export(make_predict_dat)
 export(mde)
 export(mt)
+export(plot_smooth_term)
 export(sim_Bai)
 export(summarize_ancillary_parameters)
 export(tune.bgam)
diff --git a/R/plot_smooth_term.R b/R/plot_smooth_term.R
index eeb91e9..41d9e10 100644
--- a/R/plot_smooth_term.R
+++ b/R/plot_smooth_term.R
@@ -1,29 +1,83 @@
-# Terms, a vector of charaters containing the variable names
-# plot, a logic variable, indicating if ggplots are produced
-plot_smooth_term <- function(mdl, terms, plot = TRUE){
+#' Plot non-linear functions for BHAM objects
+#'
+#' @param mdl
+#' @param terms a vector of characters containing the variable names
+#' @param plot a logic variable, indicating if ggplots are produced
+#'
+#' @return
+#' @export
+#'
+#' @examples
+#' library(glmnet); data("QuickStartExample")  # Load example data
+#' x <- QuickStartExample$x; colnames(x) <- paste0("X", 1:ncol(x))
+#' y <- QuickStartExample$y
+#'
+#' dat <- data.frame(x, y)
+#'
+#' spl_df <- data.frame(
+#'    Var = colnames(x),
+#'    Func = "s",
+#'   Args ="bs='cr', k=7"
+#' )
+#'
+#' train_sm_dat <- construct_smooth_data(spl_df, dat)
+#' train_smooth <- train_sm_dat$Smooth
+#' train_smooth_data <- train_sm_dat$data
+#' # Not implmented yet
+#' # mdl <- bgam(y~.-y,
+#' #             data = data.frame(train_smooth_data, y = y), family = "gaussian",
+#' #             prior = mde(), group = make_group(names(train_smooth_data))
+#' #             )
+#' mdl <-bamlasso(x = train_smooth_data, y = y, family = "gaussian",
+#'                group = make_group(names(train_smooth_data)))
+#' plot_smooth_term(mdl, "X3", train_smooth, min = min(x[,"X3"])-0.1, max = max(x[,"X3"]) + 0.1)
+#'
+
+plot_smooth_term <- function(mdl, term, smooth, min, max,  plot = TRUE){
+
+
+  # TODO: Implement the plotting function for other models
+  if(!("bmlasso" %in% class(mdl)))
+    stop("Not Implemented for bgam and bacoxph yet")
+
+  # TODO: check if term is included in the model, and smooth
 
   # Loop Start here
-  for(var in terms)
-  {
+  # for(var in terms)
+  # {
     # retrieve the smoothing object
-
-    # Retrive the min and max value of the variable
+    sm <- `[[`(smooth, term)
+    # Retrive the min and max value of the variable, see(min, and max argument of the funciton)
 
     # Construct the vector of possible value of the variable
+    .dat <- data.frame(seq(min, max, length.out = 200))
+    colnames(.dat) <- term
 
     # Construct the data matrix
+    ret <- mgcv::PredictMat(sm, data = .dat)
+    colnames(ret) <- create_smooth_name(sm)
 
     # Retrieve the coefficients
+    .coef <- mdl$coefficients[colnames(ret)]
 
     # Calculate the linear predictors
+    # browser()
+    .dat$value <- `[[`(.dat, term)
+    .dat$lp <- ret %*% .coef
 
     # Decide if we are plotting at the linear predictor scale or response scale
+    # TODO: implement later
 
     # Make plots
     if(plot){
       # Make ggplot
+      ggplot2::ggplot(.dat) +
+        ggplot2::geom_smooth(ggplot2::aes(x = value, y = lp))+
+      ggplot2::geom_point(ggplot2::aes(x = value, y = lp), alpha = 0.5)
+
     } else {
       # return Variable value, linear predicrtors, and response scale
+      stop("Not implemented")
     }
-  } # End Loop
+  # } # End Loop
 }
diff --git a/man/plot_smooth_term.Rd b/man/plot_smooth_term.Rd
new file mode 100644
index 0000000..2c5d592
--- /dev/null
+++ b/man/plot_smooth_term.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_smooth_term.R
+\name{plot_smooth_term}
+\alias{plot_smooth_term}
+\title{Plot non-linear functions for BHAM objects}
+\usage{
+plot_smooth_term(mdl, term, smooth, min, max, plot = TRUE)
+}
+\arguments{
+\item{plot}{a logic variable, indicating if ggplots are produced}
+
+\item{terms}{a vector of characters containing the variable names}
+}
+\value{
+
+}
+\description{
+Plot non-linear functions for BHAM objects
+}
+\examples{
+library(glmnet); data("QuickStartExample")  # Load example data
+x <- QuickStartExample$x; colnames(x) <- paste0("X", 1:ncol(x))
+y <- QuickStartExample$y
+
+dat <- data.frame(x, y)
+
+spl_df <- data.frame(
+   Var = colnames(x),
+   Func = "s",
+  Args ="bs='cr', k=7"
+)
+
+train_sm_dat <- construct_smooth_data(spl_df, dat)
+train_smooth <- train_sm_dat$Smooth
+train_smooth_data <- train_sm_dat$data
+# Not implmented yet
+# mdl <- bgam(y~.-y,
+#             data = data.frame(train_smooth_data, y = y), family = "gaussian",
+#             prior = mde(), group = make_group(names(train_smooth_data))
+#             )
+mdl <-bamlasso(x = train_smooth_data, y = y, family = "gaussian",
+               group = make_group(names(train_smooth_data)))
+plot_smooth_term(mdl, "X3", train_smooth, min = min(x[,"X3"])-0.1, max = max(x[,"X3"]) + 0.1)
+
+}
diff --git a/news.md b/news.md
index 52166a0..a56961a 100644
--- a/news.md
+++ b/news.md
@@ -1,3 +1,7 @@
+# BHAM (development version 0.2.0.9000)
+## Add non-linear plotting functions
+* Currently implemented for bamlasso models, without error prevention mechanisms.
+
 # BHAM 0.2.0
 ## Model Change
 * Add hierarchical enforcement to the bi-level selection by imposing a dependency to the non-linear components inclusion indicator prior