diff --git a/R/ModelFit.R b/R/ModelFit.R
index 385a011c..7916c4cb 100644
--- a/R/ModelFit.R
+++ b/R/ModelFit.R
@@ -223,6 +223,7 @@ fitCyclopsModel <- function(cyclopsData,
         }
 
         .cyclopsSetBeta(cyclopsData$cyclopsInterfacePtr, startingCoefficients)
+        .cyclopsSetStartingBeta(cyclopsData$cyclopsInterfacePtr, startingCoefficients)
     }
 
     if (!is.null(fixedCoefficients)) {
diff --git a/R/RcppExports.R b/R/RcppExports.R
index d4d5369b..702c104a 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -25,10 +25,22 @@
     invisible(.Call(`_Cyclops_cyclopsSetBeta`, inRcppCcdInterface, beta))
 }
 
+.cyclopsGetBeta <- function(inRcppCcdInterface, index) {
+    .Call(`_Cyclops_cyclopsGetBeta`, inRcppCcdInterface, index)
+}
+
+.cyclopsSetStartingBeta <- function(inRcppCcdInterface, inStartingBeta) {
+    invisible(.Call(`_Cyclops_cyclopsSetStartingBeta`, inRcppCcdInterface, inStartingBeta))
+}
+
 .cyclopsSetFixedBeta <- function(inRcppCcdInterface, beta, fixed) {
     invisible(.Call(`_Cyclops_cyclopsSetFixedBeta`, inRcppCcdInterface, beta, fixed))
 }
 
+.cyclopsGetFixedBeta <- function(inRcppCcdInterface, index) {
+    .Call(`_Cyclops_cyclopsGetFixedBeta`, inRcppCcdInterface, index)
+}
+
 .cyclopsGetIsRegularized <- function(inRcppCcdInterface, index) {
     .Call(`_Cyclops_cyclopsGetIsRegularized`, inRcppCcdInterface, index)
 }
diff --git a/src/RcppCyclopsInterface.cpp b/src/RcppCyclopsInterface.cpp
index ae8ff12e..c2ccab92 100644
--- a/src/RcppCyclopsInterface.cpp
+++ b/src/RcppCyclopsInterface.cpp
@@ -116,6 +116,20 @@ void cyclopsSetBeta(SEXP inRcppCcdInterface, const std::vector<double>& beta) {
     interface->getCcd().setBeta(beta);
 }
 
+// [[Rcpp::export(.cyclopsGetBeta)]]
+double cyclopsGetBeta(SEXP inRcppCcdInterface, const int index) {
+	using namespace bsccs;
+	XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+	return interface->getCcd().getBeta(index);
+}
+
+// [[Rcpp::export(.cyclopsSetStartingBeta)]]
+void cyclopsSetStartingBeta(SEXP inRcppCcdInterface, const std::vector<double>& inStartingBeta) {
+	using namespace bsccs;
+	XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+	interface->getCcd().setStartingBeta(inStartingBeta);
+}
+
 // [[Rcpp::export(.cyclopsSetFixedBeta)]]
 void cyclopsSetFixedBeta(SEXP inRcppCcdInterface, int beta, bool fixed) {
     using namespace bsccs;
@@ -124,6 +138,14 @@ void cyclopsSetFixedBeta(SEXP inRcppCcdInterface, int beta, bool fixed) {
     interface->getCcd().setFixedBeta(beta - 1, fixed);
 }
 
+// [[Rcpp::export(.cyclopsGetFixedBeta)]]
+bool cyclopsGetFixedBeta(SEXP inRcppCcdInterface, const int index){
+	using namespace bsccs;
+	XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+	return interface->getCcd().getFixedBeta(index);
+}
+
+
 // [[Rcpp::export(".cyclopsGetIsRegularized")]]
 bool cyclopsGetIsRegularized(SEXP inRcppCcdInterface, const int index) {
     using namespace bsccs;
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 67ad4bce..bf33c7c5 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -73,6 +73,29 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// cyclopsGetBeta
+double cyclopsGetBeta(SEXP inRcppCcdInterface, const int index);
+RcppExport SEXP _Cyclops_cyclopsGetBeta(SEXP inRcppCcdInterfaceSEXP, SEXP indexSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const int >::type index(indexSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetBeta(inRcppCcdInterface, index));
+    return rcpp_result_gen;
+END_RCPP
+}
+// cyclopsSetStartingBeta
+void cyclopsSetStartingBeta(SEXP inRcppCcdInterface, const std::vector<double>& inStartingBeta);
+RcppExport SEXP _Cyclops_cyclopsSetStartingBeta(SEXP inRcppCcdInterfaceSEXP, SEXP inStartingBetaSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const std::vector<double>& >::type inStartingBeta(inStartingBetaSEXP);
+    cyclopsSetStartingBeta(inRcppCcdInterface, inStartingBeta);
+    return R_NilValue;
+END_RCPP
+}
 // cyclopsSetFixedBeta
 void cyclopsSetFixedBeta(SEXP inRcppCcdInterface, int beta, bool fixed);
 RcppExport SEXP _Cyclops_cyclopsSetFixedBeta(SEXP inRcppCcdInterfaceSEXP, SEXP betaSEXP, SEXP fixedSEXP) {
@@ -85,6 +108,18 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// cyclopsGetFixedBeta
+bool cyclopsGetFixedBeta(SEXP inRcppCcdInterface, const int index);
+RcppExport SEXP _Cyclops_cyclopsGetFixedBeta(SEXP inRcppCcdInterfaceSEXP, SEXP indexSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const int >::type index(indexSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetFixedBeta(inRcppCcdInterface, index));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cyclopsGetIsRegularized
 bool cyclopsGetIsRegularized(SEXP inRcppCcdInterface, const int index);
 RcppExport SEXP _Cyclops_cyclopsGetIsRegularized(SEXP inRcppCcdInterfaceSEXP, SEXP indexSEXP) {
@@ -766,7 +801,10 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsGetUseOffsetNames", (DL_FUNC) &_Cyclops_cyclopsGetUseOffsetNames, 0},
     {"_Cyclops_cyclopsGetComputeDevice", (DL_FUNC) &_Cyclops_cyclopsGetComputeDevice, 1},
     {"_Cyclops_cyclopsSetBeta", (DL_FUNC) &_Cyclops_cyclopsSetBeta, 2},
+    {"_Cyclops_cyclopsGetBeta", (DL_FUNC) &_Cyclops_cyclopsGetBeta, 2},
+    {"_Cyclops_cyclopsSetStartingBeta", (DL_FUNC) &_Cyclops_cyclopsSetStartingBeta, 2},
     {"_Cyclops_cyclopsSetFixedBeta", (DL_FUNC) &_Cyclops_cyclopsSetFixedBeta, 3},
+    {"_Cyclops_cyclopsGetFixedBeta", (DL_FUNC) &_Cyclops_cyclopsGetFixedBeta, 2},
     {"_Cyclops_cyclopsGetIsRegularized", (DL_FUNC) &_Cyclops_cyclopsGetIsRegularized, 2},
     {"_Cyclops_cyclopsSetWeights", (DL_FUNC) &_Cyclops_cyclopsSetWeights, 2},
     {"_Cyclops_cyclopsSetCensorWeights", (DL_FUNC) &_Cyclops_cyclopsSetCensorWeights, 2},
diff --git a/src/cyclops/CyclicCoordinateDescent.cpp b/src/cyclops/CyclicCoordinateDescent.cpp
index 6e54b9f2..daa04671 100644
--- a/src/cyclops/CyclicCoordinateDescent.cpp
+++ b/src/cyclops/CyclicCoordinateDescent.cpp
@@ -207,10 +207,14 @@ void CyclicCoordinateDescent::init(bool offset) {
 	hDelta.resize(J, static_cast<double>(initialBound));
 	hBeta.resize(J, static_cast<double>(0.0));
 
+	// initialize starting betas to default value 0.0
+	startingBeta.resize(J, static_cast<double>(0.0));
+	
 // 	hXBeta.resize(K, static_cast<double>(0.0));
 // 	hXBetaSave.resize(K, static_cast<double>(0.0));
 
 	fixBeta.resize(J, false);
+	
 	hWeights.resize(0);
     cWeights.resize(0); // ESK: For censor weights
 
@@ -221,6 +225,7 @@ void CyclicCoordinateDescent::init(bool offset) {
 	varianceKnown = false;
 	if (offset) {
 		hBeta[0] = static_cast<double>(1);
+		startingBeta[0] = static_cast<double>(1);
 		fixBeta[0] = true;
 		xBetaKnown = false;
 	} else {
@@ -261,7 +266,13 @@ void CyclicCoordinateDescent::resetBeta(void) {
 
     auto start = hXI.getHasOffsetCovariate() ? 1 : 0;
     for (auto j = start; j < J; j++) {
-		hBeta[j] = 0.0;
+
+    	// check if a non-zero starting beta is present
+    	if (startingBeta[j] != 0.0) {
+    		hBeta[j] = startingBeta[j];
+    	} else {
+    		hBeta[j] = 0.0;
+    	}
 	}
 	computeXBeta();
 	sufficientStatisticsKnown = false;
@@ -474,6 +485,13 @@ void CyclicCoordinateDescent::setBeta(const std::vector<double>& beta) {
 	varianceKnown = false;
 }
 
+void CyclicCoordinateDescent::setStartingBeta(const std::vector<double>& inStartingBeta) {
+	// ToDo: This functionality could be merged into setBeta()
+	for (int j = 0; j < J; ++j) {
+		startingBeta[j] = inStartingBeta[j];
+	}
+}
+
 void CyclicCoordinateDescent::setBeta(int i, double beta) {
 #define PROCESS_IN_MS
 #ifdef PROCESS_IN_MS
diff --git a/src/cyclops/CyclicCoordinateDescent.h b/src/cyclops/CyclicCoordinateDescent.h
index 0228a192..9ccba618 100644
--- a/src/cyclops/CyclicCoordinateDescent.h
+++ b/src/cyclops/CyclicCoordinateDescent.h
@@ -133,6 +133,8 @@ class CyclicCoordinateDescent {
 //	template <typename T>
 	void setBeta(const std::vector<double>& beta);
 
+	void setStartingBeta(const std::vector<double>& inStartingBeta);
+	
 	void setBeta(int i, double beta);
 
 //	void double getHessianComponent(int i, int j);
@@ -330,6 +332,7 @@ class CyclicCoordinateDescent {
 
 	typedef std::vector<double> DoubleVector;
 	DoubleVector hBeta;
+	DoubleVector startingBeta;
 
 // 	DoubleVector& hXBeta; // TODO Delegate to ModelSpecifics
 // 	DoubleVector& hXBetaSave; // Delegate