Rho-square calculation

[Amir Ghorbani]

Hello,

I've been exploring the Biogeme source code to understand how rho-squared, null likelihood and estimated model likelihood are computed.

Could you assist me in locating this part of the code?

Best regards,
Amir

[Michel Bierlaire]

In the documentation, it is defined in Section 7 of https://transp-or.epfl.ch/documents/technicalReports/Bier23.pdf

In the code, the output of the estimation is transmitted to the “results” object.
https://github.com/michelbierlaire/biogeme/blob/master/src/biogeme/biogeme.py, rows1512-1516

The calculation of the null log likelihood is made at row 902.

The calculation of the rho bar squared is made in the file
https://github.com/michelbierlaire/biogeme/blob/master/src/biogeme/results.py
in the function “calculate_stats”, row 408.

> --
> You received this message because you are subscribed to the Google Groups "Biogeme" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to biogeme+u...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/biogeme/a999ec4d-9e34-46c1-8646-d64fcb83bf89n%40googlegroups.com.

Michel Bierlaire
Transport and Mobility Laboratory
School of Architecture, Civil and Environmental Engineering
EPFL - Ecole Polytechnique Fédérale de Lausanne
http://transp-or.epfl.ch
http://people.epfl.ch/michel.bierlaire

[Michel Bierlaire]

The likelihood itself is calculated using the C++ code.
It is distributed separately, as it requires a c++ compiler to be installed.
https://github.com/michelbierlaire/cythonbiogeme

The main function is available in this file:
https://github.com/michelbierlaire/cythonbiogeme/blob/main/src/cythonbiogeme/cpp/biogeme.cc



> On 5 Sep 2023, at 09:45, Amir Ghorbani <ghor...@student.unimelb.edu.au> wrote:
>
>
> Hi
>
> Thanks for the response.
>
> Whats the exact formula in the code for computing likelihood ? I can not see that . Is there any division by a scale ?
>
> Best regards, Amir
>
>> On 5 Sep 2023, at 3:43 pm, Amir Ghorbani <saken...@gmail.com> wrote:
>>
>> 
>>
>> Sent from my iPhone
>>
>> Begin forwarded message:
>>
>>> From: Michel Bierlaire <michel.b...@epfl.ch>
>>> Date: 5 September 2023 at 3:16:48 pm AEST
>>> To: saken...@gmail.com
>>> Cc: Michel Bierlaire <michel.b...@epfl.ch>, Biogeme <bio...@googlegroups.com>
>>> Subject: Re: [biogeme] Rho-square calculation
>>>
>>> In the documentation, it is defined in Section 7 of https://transp-or.epfl.ch/documents/technicalReports/Bier23.pdf

[Amir Ghorbani]

Still need help finding the formula. I see no logarithm in the code, which is necessary for loglikelihood calculation.

Can you please copy and paste the part of the code here?

Best regards, Amir

[Michel Bierlaire]

The formula provided by the user is the contribution to the log likelihood of each observation.

If you look at the examples (e.g. https://github.com/michelbierlaire/biogeme/blob/master/examples/swissmetro/b01logit.py), the formula that is built is the *log* of the choice model (see row 48).

> To view this discussion on the web visit https://groups.google.com/d/msgid/biogeme/ec351614-1470-429e-bba5-57bedfdd0013n%40googlegroups.com.

[Amir Ghorbani]

still can not see the formulas for constructing the likelihood function clearly. I can see nested references only. 

By the way , I'm trying to calculate the rho-squared for the baseline estimated by Biogeme but get a negative value( please see below). What could be the reason? Is there any scaling involved in biogeme ?


##Baseline (biogeme) -swiss
print('data[0]',data[0],'data_crop[0]',data_crop[0])
def calculate_probability_logit1(V_matrix):
    exp_V = np.exp(V_matrix)
    sum_exp_V = np.sum(exp_V, axis=1, keepdims=True)
    return exp_V / sum_exp_V

def log_likelihood1(y, P):
    threshold = 1e-10
    log_P = np.where(P > threshold, np.log(P), np.log(threshold))
    return np.sum(log_P * y)

def log_likelihood_null1(Ie):
    """Calculate the log-likelihood of the null model."""
    epsilon = 1e-1  # Small constant to avoid taking the log of zero
    num_alternatives = Ie.shape[1]  # Assume that Ie is an array where each row is an individual and each column is an alternative
    P_null = 1 / num_alternatives  # Equal probability for all alternatives
    print('P_null' ,P_null)
    return np.sum(np.log(P_null) * Ie)


# Variables - these are matrices
car_traveltime = data[:,25]
train_traveltime = data[:,18]
swissmetro_traveltime = data[:,21]
train_cost=data[:,19]
car_cost=data[:,26]
swissmetro_cost=data[:,22]
# Coefficients from the Swissmetro model
asc_car = -0.155
asc_train = -0.701
asc_swissmetro = 0  # Base category
b_time = -1.28
b_cost = -1.08
##scale
s=1e2
train_traveltime=train_traveltime/s
train_cost=train_cost/s
swissmetro_traveltime=swissmetro_traveltime/s
swissmetro_cost=swissmetro_cost/s
car_traveltime=car_traveltime/s
car_cost=car_cost/s
# Utility functions
V_car = asc_car + np.multiply(b_time, car_traveltime) + np.multiply(b_cost, car_cost)
V_train = asc_train + np.multiply(b_time, train_traveltime)+np.multiply(b_cost, train_cost)
V_swissmetro = asc_swissmetro + np.multiply(b_time, swissmetro_traveltime) +np.multiply(b_cost, swissmetro_cost)

V_matrix = np.column_stack((V_train,V_swissmetro, V_car))

print(V_matrix)

# Initialize a zero array of the same shape as V
V_transformed = np.zeros_like(V_matrix)

# Find the indices of the maximum values along axis 1 (columns)
max_indices = np.argmax(V_matrix, axis=1)

# Assign 1 to the maximum values in each row
for row_index, max_index in enumerate(max_indices):
    V_transformed[row_index, max_index] = 1
print(V_transformed)

Ie_pred=V_transformed
probs = calculate_probability_logit1(V_matrix)
loglik = log_likelihood1(Ie, probs)
loglik_null = log_likelihood_null1(Ie)

# Calculate rho-squared
rho_squared = 1 - (abs(loglik) / abs(loglik_null))

print('probs',probs,'loglik',loglik,'loglik_null ' ,loglik_null , 'rho_squared' , rho_squared )
num_params=16
num_observations = len(Ie)
aic, bic = calculate_aic_bic(loglik, num_params, num_observations)
print('Full data AIC for logit: ', aic)
print('Full data BIC for logit: ', bic)


# Calculate the predicted exposures using the best solution
x_opt = mean_eps_value
Ie_Tot_pred = Ie_pred
#Ie_Train_pred=np.round(sigmoid(np.dot(lib_sample, x_opt.reshape((-1, num_modes))) + epsi_draw_sample ))


#R_sq = 1 - ss_res / ss_tot
#R_sq_Train= 1 - (np.abs(Ie_sample -Ie_Train_pred)).sum() / (np.abs(Ie_sample - np.mean(Ie_sample))).sum()
R_sq_Total= 1 - (np.abs(Ie -Ie_Tot_pred )).sum() / (np.abs(Ie - np.mean(Ie))).sum()


print('R_sq_Total_logit',R_sq_Total)
RMSE_Total_logit = np.sqrt(np.mean((Ie - Ie_Tot_pred)**2))
print('RMSE_Total_logit', RMSE_Total_logit)


output:

data[0] [ 2 0 1 1 1 0 1 1 0 3 0 2 0 2 1 1 1 1 112 0 120 63 0 20 0 117 65 2] data_crop[0] [112 0 120 63 0 20 0 117 65 2] [[-2.1346 -0.8064 -2.3546] [-2.0194 -0.768 -2.5598] [-2.365 -0.8576 -2.2142] ... [-2.0834 -0.64 -1.8702] [-2.3394 -0.6784 -2.3022] [-2.0834 -0.6784 -2.299 ]] [[0. 1. 0.] [0. 1. 0.] [0. 1. 0.] ... [0. 1. 0.] [0. 1. 0.] [0. 1. 0.]] P_null 0.3333333333333333 probs [[0.17931551 0.67678042 0.14390407] [0.1969377 0.6883431 0.1147192 ] [0.14975125 0.67612338 0.17412537] ... [0.15449497 0.65429709 0.19120794] [0.1369399 0.72093009 0.14213001] [0.17002221 0.69293008 0.13704771]] loglik -13169.342032565199 loglik_null -6159.919102562092 rho_squared -1.1379082766018276 Full data AIC for logit: 26370.684065130397 Full data BIC for logit: 26476.792402668376 R_sq_Total_logit 0.3600856072766184 RMSE_Total_logit 0.5332976624418555