Najpierw zdefiniujmy nasze zmienne niezależne i zależne:

target = "revenue"
media_channels = ["tv_S", "ooh_S", "print_S", "facebook_S", "search_S"]
organic_channels = ["newsletter"]
control_features = ["trend", "season", "holiday", "competitor_sales_B", "events"]
features = control_features + media_channels + organic_channels

W sumie mamy 11 zmiennych niezależnych, 5 z nich to kanały wydatków na media plus 1 kanał organiczny, a 5 to zmienne kontrolne. Zastosujmy klasyczną wielokrotną regresję liniową przy użyciu pakietu statsmodels, a następnie sprawdźmy, czy krzywa pasuje do niej krzywa_dopasowanie daje takie same wyniki dla nieograniczonej regresji liniowej

from statsmodels.regression import linear_model 
from statsmodels import toolsdata_with_intercept = tools.add_constant(final_data[features])
ols = linear_model.OLS(endog = data[target], 
exog = data_with_intercept)
ols_res = ols.fit()
print(ols_res.summary())
print(ols_res.params)

Dokumentacja ww krzywa_dopasowanie funkcja mówi, że oczekuje funkcji modelowej:

Musi przyjąć zmienną niezależną jako pierwszy argument i parametry pasujące jako osobne pozostałe argumenty

Funkcja modelu opisuje naszą zależność liniową jako:

Funkcja modelowa powinna przyjmować tyle parametrów, ile zdefiniujemy i nie ma możliwości dostosowania tej funkcji do innej liczby parametrów. Ponieważ mamy 11 współczynników plus 1 wyraz wolny, musimy podać 12 parametrów:

def linear_function(data, a, b1, b2, b3, b4, 
b5, b6, b7, b8, b9, b10, b11):
return a + b1 * data[:, 0]\
+ b2 * data[:, 1]\
+ b3 * data[:, 2]\
+ b4 * data[:, 3]\
+ b5 * data[:, 4]\
+ b6 * data[:, 5]\
+ b7 * data[:, 6]\
+ b8 * data[:, 7]\
+ b9 * data[:, 8]\
+ b10 * data[:, 9]\
+ b11 * data[:, 10]

Dzwonić krzywa_dopasowanie dostarczając funkcję modelu, zmienne niezależne jako macierz NumPy, zmienną docelową jako tablicę NumPy i jawnie ustawiając metodę na „lm”, która wykorzystuje liniową optymalizację metodą najmniejszych kwadratów

curve_fit_coefs, _ = curve_fit(f = linear_function, 
method = "lm", 
xdata=final_data[features].values, 
ydata = final_data[target].values)#print coefficients
pd.DataFrame(curve_fit_coefs, index=["const"] + features, columns=["coefficient"])

Widzimy, że w obu przypadkach otrzymujemy te same współczynniki. Punkt przecięcia (const) i szukaj_S kanał medialny są negatywne. Ograniczmy przechwytywanie i wszystkie kanały medialne tak, aby były nieujemne. Na Githubie Robyn toczy się interesująca dyskusja na temat tego, czy przecięcia powinny być nieujemne, czy nieograniczone. Sprawdź ten wątek.

Najpierw musimy zdefiniować funkcję pomocniczą, aby szybko wygenerować ograniczenia w formacie krzywa_dopasowanie akceptuje: tablica o długości równej liczbie parametrów

miedza: 2-krotka typu array_like, opcjonalnie: Dolne i górne granice parametrów. Domyślnie brak ograniczeń. Każdy element krotki musi być albo an tablica o długości równej liczbie parametrówlub skalar (w tym przypadku przyjmuje się, że granica jest taka sama dla wszystkich parametrów). Używać np.inf z odpowiednim znakiem, aby wyłączyć ograniczenia wszystkich lub niektórych parametrów.

def prepare_bounds(intercept_value, 
control_length, 
media_length, 
control_value, 
media_value):
lower_bounds_array = []
lower_bounds_array.append(intercept_value)    for i in range(control_length):
lower_bounds_array.append(control_value)    for i in range(media_length):
lower_bounds_array.append(media_value)
return lower_bounds_array

Przygotujmy ograniczenia:

• Dolne granice dla kanałów przechwytujących i medialnych wynoszą 0 i -nieskończoność dla zmiennych kontrolnych.
• Górne granice to nieskończoność dla wszystkich zmiennych.

lower_bounds_array = \
prepare_bounds(intercept_value = 0, 
control_length = len(control_features), 
media_length = len(media_channels) + len(organic_channels), 
control_value = -np.inf, 
media_value = 0)
upper_bounds_array = \
prepare_bounds(intercept_value = np.inf, 
control_length = len(control_features), 
media_length = len(media_channels) + len(organic_channels), 
control_value = np.inf, 
media_value = np.inf)

Dolne granice to:

lower_bounds_array
#[0, -inf, -inf, -inf, -inf, -inf, 0, 0, 0, 0, 0, 0]

Górne granice to:

upper_bounds_array
#[inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf, inf]

Aplikujmy krzywa_dopasowanie z ograniczeniami

curve_fit_bounded_coefs, _ = \
curve_fit(f = linear_function,
xdata=final_data[features].values, 
ydata = final_data[target].values,
bounds = (tuple(lower_bounds_array), 
tuple(upper_bounds_array)))#print coefficients
pd.DataFrame(curve_fit_bounded_coefs, index = ["const"] + features, columns = ["coefficient"])

Tym razem wszystkie parametry modelu są dodatnie.

To rozwiązanie ma kilka wad:

• Funkcja modelu wymaga wcześniejszego określenia dokładnej liczby parametrów, co utrudnia to w scenariuszach, gdy chcemy poeksperymentować z inną liczbą parametrów
• Klasyczna regresja liniowa jest podatna na nadmierne dopasowanie, nie działa dobrze ze zmiennymi skorelowanymi lub nadmiarowymi i nie pozwala dobrze uogólniać. W modelowaniu marketingu mix często zdarza się, że zmienne są skorelowane lub gdy dane są ograniczone i wymagane są lepsze podejścia do modelowania, takie jak regresja Ridge’a. Regresję grzbietową zastosowano w rzeczywistych scenariuszach, aby przezwyciężyć wyzwania związane z nadmiernym dopasowaniem i współliniowością za pomocą techniki regularyzacji, która karze współczynniki poprzez zmniejszenie ich wielkości. Stopień penalizacji kontrolowany jest za pomocą parametru lambda.

Dlaczego zdecydowałem się użyć R glmnet w Pythonie? Po prostu nie mogłem znaleźć lepszego rozwiązania. Mógłbym znaleźć opakowania glmnet w Pythonie, ale ze względu na zależności w Fortranie nie mogłem ich skompilować na moim komputerze z systemem Windows. Poza tym glmnet ma od razu dwie przydatne funkcje: cv.glmnet która przeprowadza weryfikację krzyżową i określa optymalny parametr lambda, oraz funkcję glmnet, która buduje ostateczny model. Obie funkcje dokonują standaryzacji danych i pozwalają kontrolować znak współczynników i wyraz wolny. Musiałem tylko wymyślić, jak wywołać kod R z Pythona. Na szczęście istnieje pakiet RPy2, interfejs Pythona do języka R, przeznaczony właśnie do tego:

rpy2 uruchamia osadzony język R, zapewniając dostęp do niego z Pythona przy użyciu własnego C-API języka R poprzez:

interfejs wysokiego poziomu, dzięki któremu funkcje R są obiektami podobnymi do funkcji Pythona i zapewniają płynną konwersję do struktur danych numpy i pandas

interfejs niskiego poziomu bliższy C-API

Instalowanie R i RPy2

Oczywiście powinieneś mieć zainstalowany R. Następnie zainstaluj RPy2:

pip install rpy2

Sprawdź, czy RPy2 rozpoznaje Twoją instalację R, importując RPy2:

from rpy2 import robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
import rpy2.situation
base = importr("base")
glm = importr("glmnet")

Jeśli import przebiegł pomyślnie, uzyskaj informacje na temat instalacji języka R:

for row in rpy2.situation.iter_info():
print(row)

Jeśli podczas importowania wystąpiły błędy, możesz spróbować wykonać następujące kroki:

• sprawdź, czy itp/Rprofil plik nie ładuje żadnej biblioteki
• Stwórz R_DOM środowisku ze ścieżką do instalacji R

Zawijanie pakietu glmnet

Proces owijania składa się z trzech etapów

• Napisanie funkcji w R, która wywołuje regresję Ridge’a i zwraca interesujący nas wynik
• Pisanie konwerterów typu Python na R dla Pandas DataFrames, Boolean i Float
• Łączenie wszystkiego w funkcji Pythona, która przygotowuje typy danych, wywołuje kod R i zwraca wyniki

Krok 1: Napisanie kodu R: dopasowanie regresji Ridge’a w R

Zapewnia RPy2 rpy2.robjects.r funkcja, która pobiera dowolny ciąg znaków i ocenia go jako kod R:

run_glmnet_in_r = \
ro.r("""
#pure R code
function(x_train, 
y_train, 
x_test, 
lambda, 
is_intercept, 
lower_limits, 
upper_limits, 
alpha = 0){ 
mod = glmnet(x_train,
y_train,
family = 'gaussian',
alpha = alpha,
lambda = lambda, 
intercept = is_intercept,
lower.limits = lower_limits,
upper.limits = upper_limits,
type_measure = 'mse')coefs = as.data.frame(as.matrix(coef(mod)))
names(coefs)[1] = "value"
y_pred = predict(mod, s = lambda, newx = x_test)
return (list(coefs = coefs, y_pred = y_pred))
}""")

Zdefiniowałem czystą funkcję R, która pasuje do regresji Ridge’a i oczekuje 8 parametrów:

• x_train — zestaw treningowy
• y_train — odpowiedź/zmienna docelowa
• x_test — zbiór testowy, którego przewidywania otrzymujemy z dopasowanego modelu
• lambda — wartość lambda, z której otrzymujemy cv.glmnet
• is_intercept — czy chcemy, aby wyraz wolny był dopasowany (nieograniczony) czy ustawiony na 0
• dolne_limity, górne_limity — wektor dolnych/górnych limitów dla każdej zmiennej
• alpha — glmnet dopasowuje regresje Ridge, Lasso lub Elasticnet kontrolowane przez alfa. Gdy alfa 0 — zamontowana jest regresja grzbietu.

Zwracana jest lista współczynników i przewidywań zestawu testowego.

Krok 2: Pisanie konwerterów typu Python na R dla Pandas DataFrames, Boolean i Float

Poniższy kod konwertuje ramkę danych Pandas na ramkę danych R:

with ro.conversion.localconverter(ro.default_converter + pandas2ri.converter):
r_df = ro.conversion.py2rpy(x_train)

Jednak glmnet oczekuje macierzy jako danych wejściowych. Musiałem napisać własny konwerter, który pobiera ramkę danych R i konwertuje ją na macierz:

data_frame_to_r_matrix = ro.r('function(x) as.matrix(x)')

Wektory wartości docelowych, dolnych i górnych granic są konwertowane za pomocą RPy2:

r_y_train = ro.vectors.FloatVector(y_train)
r_lower_limits = ro.vectors.FloatVector(lower_limits)
r_upper_limits = ro.vectors.FloatVector(upper_limits)

Wartość skalarna zmiennoprzecinkowa lambda jest konwertowana tak, jakby była wektorem o pojedynczej wartości:

lambda_best_float = ro.vectors.FloatVector([lambda_best])

Zmienna logiczna jest_przechwytem konwertuje się podobnie:

is_intercept_bool = ro.vectors.BoolVector([is_intercept])

Konwersja listy R do słownika Pythona wymaga niestandardowego rozwiązania:

def convert_r_list_into_dictionary(r_list):
py_dict = {}
#for name in coefs.names:
keys = r_list.names
for i in range(len(keys)):
if isinstance(r_list[i], ro.vectors.DataFrame):
with ro.conversion.localconverter(ro.default_converter + pandas2ri.converter):
py_dict[keys[i]]=r_list[i]
elif isinstance(r_list[i], ro.vectors.FloatVector):
array = np.array(r_list[i])
if len(array) == 1:
array = array[0]
py_dict[keys[i]] = array
else:
py_dict[keys[i]] = r_list[i]
return py_dict

Krok 3: Połączenie wszystkich kroków w funkcji Pythona, która przygotowuje typy danych, wywołuje kod R i zwraca wyniki

def run_glmnet(x_train, 
y_train, 
x_test, 
lambda_best, 
is_intercept, 
lower_limits, 
upper_limits, 
alpha = 0):
#type conversions
with ro.conversion.localconverter(ro.default_converter + pandas2ri.converter):
r_df = ro.conversion.py2rpy(x_train)
r_test_df = ro.conversion.py2rpy(x_test)
r_x_train = data_frame_to_r_matrix(r_df)
r_x_test = data_frame_to_r_matrix(r_test_df)
r_y_train = ro.vectors.FloatVector(y_train)
r_lower_limits = ro.vectors.FloatVector(lower_limits)
r_upper_limits = ro.vectors.FloatVector(upper_limits)
lambda_best_float = ro.vectors.FloatVector([lambda_best])
is_intercept_bool = ro.vectors.BoolVector([is_intercept])
#call glmnet
mod = run_glmnet_in_r(r_x_train, 
r_y_train, 
r_x_test, 
lambda_best_float, 
is_intercept_bool, 
r_lower_limits, 
r_upper_limits, 
alpha)
    #prepare the results
mod = convert_r_list_into_dictionary(mod)
mod["coefs"] = mod["coefs"].reset_index()
return mod

Dopasowanie regresji kalenicy

Zdefiniujmy dolną i górną granicę:

# Set lower and upper bounds for feature coefficients
lower_limits = np.array([-np.inf for _ in range(len(control_features))] + [0 for _ in range(len(media_channels) + len(organic_channels))])
upper_limits = np.array([np.inf for _ in range(len(control_features))] + [np.inf for _ in range(len(media_channels) + len(organic_channels))])print(lower_limits)
print(upper_limits)#[-inf -inf -inf -inf -inf   0.   0.   0.   0.   0.   0.]
#[inf inf inf inf inf inf inf inf inf inf inf]

Tym razem definiujemy ograniczenia tylko dla zmiennych. Przecięcie jest jawnie kontrolowane przez parametr funkcji.

Następnie sprawdzamy krzyżowo, aby znaleźć najlepszy parametr lambda:

cv_glmnet = run_cv_glmnet(x_train = final_data[features], 
y_train = final_data[target], 
is_intercept = True, 
lower_limits = lower_limits , 
upper_limits = upper_limits)print(cv_glmnet["lambda_best"])
#166018.4312548283

Na koniec dopasowujemy regresję Ridge’a za pomocą znalezionej lambdy:

results = run_glmnet(x_train = final_data[features], 
y_train = final_data[target], 
x_test = final_data[features], 
lambda_best = cv_glmnet["lambda_best"], 
is_intercept = True, 
lower_limits = lower_limits, 
upper_limits = upper_limits)results.keys()
#dict_keys(['coefs', 'y_pred'])

Ostateczne współczynniki to:

results["coefs"]

W niektórych sytuacjach, kierując się logiką biznesową jak w Modelowaniu Marketing Mix, chcemy wymusić regresję liniową, aby znaleźć rozwiązanie o dodatnich współczynnikach. Jednak nie wszystkie pakiety statystyczne obsługują tę funkcję i potrzebne jest pewne obejście.

Pokazałem dwa sposoby ograniczania współczynników za pomocą nieliniowej funkcji optymalizacji metodą najmniejszych kwadratów SciPy krzywa_dopasowanie i regresja Ridge’a z pakietu R glmnet. To drugie podejście wymaga zawinięcia kodu R przy użyciu interfejsu RPy2 i umożliwia wywoływanie dowolnych funkcji R z Pythona.

Cały kod można pobrać z mojego repozytorium Github

Dziękuje za przeczytanie!