diff --git a/yli/regress.py b/yli/regress.py
index 4470ae0..1e281eb 100644
--- a/yli/regress.py
+++ b/yli/regress.py
@@ -102,6 +102,7 @@ def regress(
 	*,
 	nan_policy='warn',
 	exposure=None,
+	method=None, maxiter=None, start_params=None,
 	bool_baselevels=False, exp=None
 ):
 	"""
@@ -119,6 +120,9 @@ def regress(
 	:type nan_policy: str
 	:param exposure: Column in *df* for the exposure variable (numeric, some models only)
 	:type exposure: str
+	:param method: See statsmodels *model.fit*
+	:param maxiter: See statsmodels *model.fit*
+	:param start_params: See statsmodels *model.fit*
 	:param bool_baselevels: Show reference categories for boolean independent variables even if reference category is *False*
 	:type bool_baselevels: bool
 	:param exp: Report exponentiated parameters rather than raw parameters, default (*None*) is to autodetect based on *model_class*
@@ -135,11 +139,19 @@ def regress(
 	df_ref = weakref.ref(df)
 	dmatrices, dep_categories = df_to_dmatrices(df, dep, formula, nan_policy)
 	
+	# Build function call arguments
+	fit_kwargs = {}
+	if exposure is not None:
+		fit_kwargs['exposure'] = exposure
+	if method is not None:
+		fit_kwargs['method'] = method
+	if maxiter is not None:
+		fit_kwargs['maxiter'] = maxiter
+	if start_params is not None:
+		fit_kwargs['start_params'] = start_params
+	
 	# Fit model
-	if exposure is None:
-		result = model_class.fit(dmatrices[0], dmatrices[1])
-	else:
-		result = model_class.fit(dmatrices[0], dmatrices[1], exposure=exposure)
+	result = model_class.fit(dmatrices[0], dmatrices[1], **fit_kwargs)
 	
 	# Fill in general information
 	result.df = df_ref
@@ -1131,13 +1143,13 @@ class Poisson(RegressionModel):
 		return 'Poisson'
 	
 	@classmethod
-	def fit(cls, data_dep, data_ind, exposure=None):
+	def fit(cls, data_dep, data_ind, exposure=None, method='newton', maxiter=None, start_params=None):
 		result = cls()
 		result.exp = True
 		result.cov_type = 'nonrobust'
 		
 		# Perform regression
-		raw_result = sm.Poisson(endog=data_dep, exog=data_ind, exposure=exposure, missing='raise').fit(disp=False)
+		raw_result = sm.Poisson(endog=data_dep, exog=data_ind, exposure=exposure, missing='raise').fit(disp=False, method=method, start_params=start_params)
 		
 		result.dof_model = raw_result.df_model
 		result.dof_resid = raw_result.df_resid