Rama Reksotinoyo

Posted on Oct 8, 2023

Intepretasi Regresi Logistik: Dari Logit Hingga Penurunan Gradient.

#statistics #computerscience #programming

Dengan menyebut yang dipertuan agung programmer andal di luar sana dan ChatGPT, saya haturkan salam sebelum menuliskan kode program pada artikel ini.

Regresi logistik merupakan satu dari banyak algoritma machine learning yang sering dipakai untuk kasus klasifikasi. Intepretasi model regresi logistik pada curva membentuk melengkung/non-linear. Ini adalah kurva regresi logistik.

Model diatas dapat diintepretasikan sebagai garis yang memilki kemiringan untuk βπ(x)[1-π(x)] terhadap parameter β.

Jika persamaan diatas dielaborasi lebih dalam lagi, π(x) menggambarkan probabilitas suatu peristiwa terjadi pada titik x. Persamaan ini mengasumsikan bahwa probabilitas peristiwa tersebut dapat diprediksi atau diperkirakan menggunakan fungsi π(x). Misalnya, jika π(x) memiliki nilai 0,5, maka probabilitas peristiwa terjadi pada titik x adalah 0,5. Dan titik 0,5 sendiri merupakan titik yang paling curam pada kurva.

Masuk ke bagian yang paling penting, yaitu intepretasi dari model yang dihasilkan. Hasilnya akan direpresentasikan dalam bentuk logit atau log odds, dengan odds sendiri adalah peluang sukses dibagi dengan peluang gagal. Dengan persamaan dibawah ini:

log(odds) = β₀ + β₁x
Intepretasi : setiap unit pertambahan pada variabel x akan mengakibatkan pertambahan pada β1 unit pada log odds.

Selanjutnya, untuk mengintepretasikan oddds, maka akan ditransformasi nilai log dari odds menjadi odds.
Dengan demikian,

odds = e^(β₀ + β₁x) = e^(β₀) * (e^(β₁))^x

Interpretasi: odds perkiraan keberhasilan akan dikalikan dengan exp(β₁) untuk setiap peningkatan satu unit pada x.

Masuk kedalam kode program, disini saya mencoba mengimplementasikan regresi logistik dengan menggunakan bahasa pemrograman C yang tentu saja dibantu oleh ChatGPT. Dataset yang saya gunakan adalah dataset mengenai horseshoe crab atau masyarat Indonesia lebih mengenal dengan sebutan belangkas. Dataset berisi satu fitur atau prediktor berupa lebar cangkang, dan target kelas atau respon berisi punya atau tidaknya satelite atau tanduk pada kepiting ini. Jumlah record pada dataset adalah 173 baris dengan data sudah bersih dari data kosong, dublikat, dan ketidak konsistenan data.

Pertama yang saya lakukan adalah membuat struct yang nantinya akan direpresentasikan sebagai objek, baik itu sebagai objek data atau pun ndarray (N-Dimensional Array) atau array berdimensi tinggi.

typedef struct {
  unsigned int num_rows;
  unsigned int num_cols;
  float **data;
  int is_square;
} Ndarray;

Lalu saya membuat fungsi yang digunakan untuk memuat dataset yang nantinya digunakan pada program ini.

Ndarray* readCSV(const char* filename, int maxRecords)
{
  FILE* file = fopen(filename, "r");
  if (file == NULL)
  {
    printf("Error opening file.\n");
    return NULL;
  }

  int records = 0;
  int read = 0;
  float** data = malloc(maxRecords * sizeof(float*));
  if (data == NULL) {
    printf("Error allocating memory.\n");
    fclose(file);
    return NULL;
  }

  while (records < maxRecords)
  {
    data[records] = malloc(2 * sizeof(float));
    if (data[records] == NULL) {
      printf("Error allocating memory.\n");
      for (int i = 0; i < records; i++) {
        free(data[i]);
      }
      free(data);
      fclose(file);
      return NULL;
    }

    read = fscanf(file, "%f,%f\n", &data[records][0], &data[records][1]);

    if (read == 2)
      records++;
    else if (read != 2 && !feof(file))
    {
      printf("File format incorrect.\n");
      for (int i = 0; i < records; i++) {
        free(data[i]);
      }
      free(data);
      fclose(file);
      return NULL;
    }

    if (feof(file))
      break;
  }

  fclose(file);

  Ndarray* crabs = malloc(sizeof(Ndarray));
  if (crabs == NULL) {
    printf("Error allocating memory.\n");
    for (int i = 0; i < records; i++) {
      free(data[i]);
    }
    free(data);
    return NULL;
  }

  crabs->num_rows = records;
  crabs->num_cols = 2;
  crabs->data = data;

  return crabs;
}

Demi mengurangi drama type casting, saya terpaksa membuat varibel yang berupa bilangan real saya rubah ke desimal.

Dilanjutkan dengan saya membuat fungsi untuk mengetahui berapa jumlah sukses(y=1) dan jumlah gagal pada dataset(y=0).

int get_n_success(Ndarray* crabs)
{
  int n_success = 0;

  for (int i = 0; i < crabs->num_rows; i++)
  {
    if (crabs->data[i][0] == 1.0)  // Perubahan pada baris ini
    {
      n_success += 1;
    }
  }

  return n_success;
}

int get_n_failed(Ndarray* crabs)
{
  int n_failed = 0;

  for (int i = 0; i < crabs->num_rows; i++)
  {
    if (crabs->data[i][0] == 0.0)
    {
      n_failed += 1;
    }
  }

  return n_failed;
}

Hasil dari kode program diatas adalah n_success sebanyak 111, dan n_failure sebanyak 62.

Selanjutnya adalah gradient decent, yang saya yakini dengan subjektifitassaya merupakan salah satu algoritma terbaik yang pernah diciptakan manusia. Yang mana, fungsi ini akan melakukan optimisasi parameter, pada kasus ini parameter hanya lebar cangkang.

Karena tujuan dari gradient descent adalah memperbarui paramenter, maka untuk mencari nilai optimum dari β₀ and β₁, memaksimalkan fungsi likelihood/log-likelihood, atau meminimalkan negative log-likelihood:
$-log(L(β)) = -Σ{i=1}^{n}(y{i}log[π(x_{i})]+(1-y_{i})log[1-π(x_{i})])$
Dan jika parameternya sebanya j, maka:

Untuk mencari parameter βj atau parameter sebanyak j maka:

Selajutnya sebelum membangun fungsi gradient decent, perlu disiapkan juga fungsi non-linear sigmoid, untuk mengintepretasi hasil klasifikasi, entah hasilnya akan masuk ke kelas 0 atau 1. Dibawah ini adalah implementasi dari fungsi sigmoid:

double sigmoid(double *X, int num_cols, double b0, double *b1) {
    double logit = b0;
    double exp_val, pi;

    for (int i = 0; i < num_cols; i++) {
        logit += X[i] * b1[i];
    }

    exp_val = exp(logit);
    pi = exp_val / (1 + exp_val);

    return pi;
}

Dilanjutkan dengan membangun fungsi yang dapat menghitung nilai sigmoid disetiap data point:

Ndarray* calculate_sigmoid_list(Ndarray* matrix, double b0, double *b1) {
    Ndarray* sigmoid_list = malloc(sizeof(Ndarray));
    sigmoid_list->num_rows = matrix->num_rows;
    sigmoid_list->num_cols = 1;
    sigmoid_list->data = malloc(sigmoid_list->num_rows * sizeof(double*));

    for (int i = 0; i < sigmoid_list->num_rows; i++) {
        sigmoid_list->data[i] = malloc(sigmoid_list->num_cols * sizeof(double));

        double* matrix_row = (double*) matrix->data[i];

        sigmoid_list->data[i][0] = sigmoid(matrix_row, matrix->num_cols, b0, b1);
    }

    return sigmoid_list;
}

Langkah selanjutnya adalah membuat fungsi kerugian atau nilai error, pada impelemtasinya, saya menerapkan negative log loss-likelihood. Dengan persamaan:
$-log(L(β)) = -Σ{i=1}^{n}(y{i}log[π(x_{i})]+(1-y_{i})log[1-π(x_{i})])$
Dan fungsi pada kode programnya adalah:

double cost_function(Ndarray* matrix, double *pi) {
    double eps = 1e-10;
    double log_loss = 0.0;

    for (int i = 0; i < matrix->num_rows; i++) {
        double log_like_success = matrix->data[i][0] * log(pi[i] + eps);
        double log_like_failure = (1 - matrix->data[i][0]) * log(1 - pi[i] + eps);
        double log_like_total = log_like_success + log_like_failure;
        log_loss -= log_like_total;
    }

    return log_loss;
}

Hasilnya adalah saya mendapat nilai kerugian sebesar:
119.91446220227053

Lalu saya menyiapkan fungsi untuk menghitung turunan dari β₀ atau intercept. Dengan persamaan:
$grad_β₀ = (\sum_{i=1}^{n} [\pi(x_{i}) - y_{i}])𝑥1$

Kode programnya adalah:

float gradient_b0(Ndarray* matrix, double *pi){
  float grad;
  for(int i=0; i<matrix->num_rows; i++){
    grad += (pi[i] - matrix->data[i][0]);
  }

  return grad;
}

Sedangkan persamaan matematis dan fungsi untuk menghitung gradient dari β₁ adalah:
$grad_β₁ = (\sum_{i=1}^{n} [\pi(x_{i}) - y_{i}])𝑥$

double* gradient_b1(Ndarray* matrix, double* pi, double* y, double b0, double* b1) {
    int num_rows = matrix->num_rows;
    int num_cols = matrix->num_cols;
    double* grad_b1 = (double*)malloc(num_cols * sizeof(double));

    for (int j = 0; j < num_cols; j++) {
        grad_b1[j] = 0.0;
        for (int i = 0; i < num_rows; i++) {
            double* matrix_row = (double*)matrix->data[i];
            double sigmoid_val = sigmoid(matrix_row, num_cols, b0, b1);
            grad_b1[j] += matrix_row[j] * (sigmoid_val - y[i]);
        }
    }

    return grad_b1;
}

Setelah fungsi menghitung gradient dari β₀ dan β₁ dibuat, langkah terakhir dalam gradient descent adalah menerapkan langkah-langkah dibawah ini:

Algoritma Optimasi:

Definisikan variable response y and variabel prediktor X.
Inisialisasi parameter yang akan optimasi β₀=0.0 and β₁=0.0.
Hitung gradient dari ∂NLL/∂𝛽𝑗, sehingga mendapatkan gradient dari β₀ and β₁.
Perbarui paremeter β₀ and β₁ $βj^{new} = βj^{old} - lr *$ ∂NLL/∂𝛽𝑜𝑙𝑑𝑗
Ulangi langkah 2-4: Δ𝛽𝑗 < nilai yang ditoleransi, or ∇𝛽𝑗NLL < nilai yang ditoleransi.

Kode program dari langkah-langka diatas adalah:

void gradient_descent(double *X, double *y, int num_rows, int num_cols, double eta, double tol, double *b0, double *b1, int *iterations, double *log_loss_list) {
    srand(time(NULL));
    double b1_initial[2] = {((double)rand() / RAND_MAX) * 0.01, ((double)rand() / RAND_MAX) * 0.01};

    // Make a criteria to run the iteration
    int continue_iteration = 1;

    // Running the iteration
    int i = 0;
    while (continue_iteration) {
        // Update i
        i++;

        // Calculate success probability (pi) from the current b0 and b1
        Ndarray matrix;
        matrix.num_rows = num_rows;
        matrix.num_cols = num_cols;
        matrix.data = malloc(matrix.num_rows * sizeof(double*));
        for (int j = 0; j < matrix.num_rows; j++) {
            matrix.data[j] = malloc(matrix.num_cols * sizeof(double));
            for (int k = 0; k < matrix.num_cols; k++) {
                matrix.data[j][k] = X[j * matrix.num_cols + k];
            }
        }
        Ndarray* sigmoid_list = calculate_sigmoid_list(&matrix, *b0, b1_initial);
        double *pi = malloc(num_rows * sizeof(double));
        for (int j = 0; j < num_rows; j++) {
            pi[j] = sigmoid_list->data[j][0];
        }

        // Calculate log loss from the current b0 and b1
        double log_loss = cost_function(&matrix, pi);
        log_loss_list[i - 1] = log_loss;

        // Calculate gradient of b0 and b1
        double grad_b0 = gradient_b0(&matrix, pi);
        double* grad_b1 = gradient_b1(&matrix, pi, y, *b0, b1_initial);

        // Update b0 and b1
        *b0 -= eta * grad_b0;
        for (int j = 0; j < num_cols; j++) {
            b1[j] -= eta * grad_b1[j];
        }

        // Check convergence
        if (i > 1 && fabs(log_loss_list[i - 1] - log_loss_list[i - 2]) < tol) {
            continue_iteration = 0;
        }

        // Free memory
        free(pi);
        free(grad_b1);
        free(sigmoid_list->data);
        free(sigmoid_list);
        for (int j = 0; j < matrix.num_rows; j++) {
            free(matrix.data[j]);
        }
        free(matrix.data);
    }

    // Store the number of iterations
    *iterations = i;

    // Print final b0 and b1
    printf("Nilai akhir b0: %f\n", *b0);
    for (int j = 0; j < num_cols; j++) {
        printf("Nilai akhir b1[%d]: %f\n", j, b1[j]);
    }
}

Hasil dari nilai optimum dari β₀ dan β₁ berturut-turut adalah:
-7.971823 dan 0.009542.

Sehingga logit model dengan parameter yang optimum adalah:
log(odds(𝜋(𝑥))=-7.971823+0.009542(𝑥)

Kesimpulan dari model logit diatas adalah:

Odds Ratio: Koefisien regresi logistik (0.009542 dalam persamaan ini) mengukur perubahan dalam log-odds dari variabel prediktor (𝑥) terhadap kemungkinan kejadian sukses (atau kategori positif) dalam masalah klasifikasi biner. Dalam hal ini, setiap peningkatan satu unit dalam 𝑥 akan menghasilkan peningkatan sebesar 0.009542 dalam log-odds dari kejadian sukses. Dalam kasus Anda, variabel 𝑥 mungkin adalah salah satu fitur atau prediktor yang digunakan dalam model.
Intercept (Konstanta): Intercept (-7.971823 dalam persamaan ini) adalah nilai log-odds ketika semua variabel prediktor lainnya adalah nol. Dalam interpretasi praktis, ini mewakili dasar atau tingkat dasar log-odds kejadian sukses ketika tidak ada pengaruh dari variabel prediktor lainnya.
Log-Odds: Persamaan ini menggambarkan logaritma dari odds (rasio peluang) bahwa kejadian sukses terjadi. Dalam beberapa kasus, Anda mungkin perlu mengonversi log-odds menjadi peluang atau probabilitas jika Anda ingin mendapatkan interpretasi yang lebih intuitif.
Pengaruh Prediktor: Koefisien regresi (0.009542) mengindikasikan seberapa kuat pengaruh 𝑥 terhadap peluang kejadian sukses. Jika nilai koefisien positif, maka peningkatan dalam 𝑥 akan meningkatkan log-odds kejadian sukses, sementara jika nilai koefisien negatif, maka peningkatan dalam 𝑥 akan mengurangi log-odds kejadian sukses.

Dapat dihitung nilai odds ratio e^0.009542 untuk mengetahui seberapa besar efek 𝑥 terhadap peluang kejadian sukses. Jika odds ratio lebih besar dari 1, maka peningkatan dalam 𝑥 akan meningkatkan peluang kejadian sukses. Jika odds ratio kurang dari 1, maka peningkatan dalam 𝑥 akan mengurangi peluang kejadian sukses.

Kode program lengkap saya tulis pada github

Saran pengembangan selanjutnya:

Visualisasi hasil gradient descent.
Menghitung interval probability.
Melakukan significance testing terhadap kemungkinan seekor kepiting memiliki satelit tidak bergantung pada lebarnya atau predictornya.
Melakukan standarisasi pada fitur atau predictor.

Referensi:
Categorical Statistic Course Pacmann

Writing your own linear algebra matrix library in C

An Introduction to Statistical Learning