a & b. Input Data dan Transformasi \(x_i\)

Berdasarkan nilai \(k\) dari 4 digit terakhir NIM:\[\begin{aligned} k &= 1 + 0 + 9 + 5 = 15 \\ 0.05 \times k &= 0.05 \times 15 = 0.75 \end{aligned} \]Setiap nilai \(x_i\) ditambahkan dengan \(0.75\). Nilai awal (inisialisasi) penduga parameter adalah: \[\beta_0^{(0)} = 1.15 \quad \text{dan} \quad \beta_1^{(0)} = 2.15\]

k <- 15

# Data Awal
xi <- c(0.5, 0.4, 0.9, 0.3, 0.6, 0.2)
ni <- c(10, 8, 12, 9, 14, 13)
yi <- c(7, 6, 8, 6, 9, 10)

# Transformasi x
x <- xi + (0.05 * k)

# Tampilan Data
df_data <- data.frame(
  "i"        = 1:length(xi),
  "xi_awal"  = xi,
  "xi_baru"  = x,
  "ni"       = ni,
  "yi"       = yi,
  "proporsi" = round(yi / ni, 4)
)

knitr::kable(df_data, 
             col.names = c("i", "$x_i$ Awal", "$x_i$ Baru ($x_i + 0.75$)", "$n_i$", "$y_i$", "Proporsi ($y_i/n_i$)"),
             align = "c", 
             caption = "Tabel Transformasi Data")

c&d. Iterasi Fisher-Scoring

Dalam regresi logistik binomial, pembaruan parameter \(\boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}\) menggunakan algoritma Fisher-Scoring dirumuskan sebagai berikut:

\[\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + \mathbf{I}^{-1}\left(\boldsymbol{\beta}^{(t)}\right) \mathbf{U}\left(\boldsymbol{\beta}^{(t)}\right)\]

Keterangan Komponen:

• Matriks Desain (\(\mathbf{X}\)): Matriks berukuran \(n \times 2\) dengan kolom pertama bernilai \(1\) (intersep) dan kolom kedua bernilai \(x_i\).
• Peluang Sukses (\(\pi_i\)): Probabilitas dugaan fungsi logit: \[\pi_i = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_i)}}\]
• Matriks Bobot (\(\mathbf{W}\)): Matriks diagonal berukuran \(n \times n\): \[\mathbf{W} = \text{diag}\left(n_i \pi_i (1 - \pi_i)\right)\]
• Matriks Informasi Fisher (\(\mathbf{I}(\boldsymbol{\beta})\)): \[\mathbf{I}(\boldsymbol{\beta}) = \mathbf{X}^T \mathbf{W} \mathbf{X}\]
• Vektor Skor (\(\mathbf{U}(\boldsymbol{\beta})\)): \[\mathbf{U}(\boldsymbol{\beta}) = \mathbf{X}^T (\mathbf{y} - \mathbf{n}\boldsymbol{\pi})\]
• Galat Baku (Standard Error): \[SE(\hat{\beta}_j) = \sqrt{\left[\mathbf{I}^{-1}\right]_{jj}}\]

# Matriks Desain X dengan Intersep
X <- cbind(1, x)
colnames(X) <- c("Intercept", "X")

# Nilai Awal Parameter (beta0 = 1.15, beta1 = 2.15)
beta <- c(1.15, 2.15)

# Fungsi Log-Likelihood
loglik <- function(b) {
  eta <- as.vector(X %*% b)
  sum(yi * eta - ni * log1p(exp(eta)))
}

jumlah_iterasi <- 5
hasil_matriks  <- list()
hasil_se       <- list()
hasil_beta     <- list()
hasil_U        <- list()

for (iter in 1:jumlah_iterasi) {
  eta <- as.vector(X %*% beta)
  p   <- plogis(eta)
  
  w <- ni * p * (1 - p)
  W <- diag(w)
  I <- t(X) %*% W %*% X
  
  U     <- t(X) %*% (yi - ni * p)
  I_inv <- solve(I)
  SE    <- sqrt(diag(I_inv))
  
  delta <- as.vector(I_inv %*% U)
  
  langkah   <- 1
  beta_baru <- beta + langkah * delta
  
  while (loglik(beta_baru) < loglik(beta) && langkah > 1e-8) {
    langkah   <- langkah / 2
    beta_baru <- beta + langkah * delta
  }
  
  hasil_matriks[[iter]] <- I
  hasil_se[[iter]]      <- SE
  hasil_beta[[iter]]    <- beta_baru
  hasil_U[[iter]]       <- U
  
  beta <- beta_baru
}

# Menyusun tabel ringkasan hasil 5 iterasi
tabel_iterasi <- data.frame(
  "Iterasi"       = 1:jumlah_iterasi,
  "Beta_0"        = round(sapply(hasil_beta, `[`, 1), 6),
  "Beta_1"        = round(sapply(hasil_beta, `[`, 2), 6),
  "SE_Beta0"      = round(sapply(hasil_se, `[`, 1), 6),
  "SE_Beta1"      = round(sapply(hasil_se, `[`, 2), 6)
)

knitr::kable(
  tabel_iterasi,
  col.names = c("Iterasi", "$\\hat{\\beta}_0$", "$\\hat{\\beta}_1$", "$SE(\\hat{\\beta}_0)$", "$SE(\\hat{\\beta}_1)$"),
  align = "c",
  caption = "Ringkasan Hasil Iterasi Fisher-Scoring"
)

e. Makna Setiap Elemen Matriks Information Fisher

Matriks Informasi Fisher \(\mathbf{I}(\boldsymbol{\beta})\) berukuran \(2 \times 2\) memiliki bentuk matematis:\[\mathbf{I}(\boldsymbol{\beta}) = \mathbf{X}^T \mathbf{W} \mathbf{X} = \begin{bmatrix} \sum_{i=1}^n n_i \pi_i (1-\pi_i) & \sum_{i=1}^n n_i \pi_i (1-\pi_i) x_i \\ \sum_{i=1}^n n_i \pi_i (1-\pi_i) x_i & \sum_{i=1}^n n_i \pi_i (1-\pi_i) x_i^2 \end{bmatrix} = \begin{bmatrix} I_{11} & I_{12} \\ I_{21} & I_{22} \end{bmatrix}\] Penjelasan tiap elemennya:

• Elemen \(I_{11} = \sum n_i \pi_i (1-\pi_i)\): Menunjukkan total jumlah informasi (bobot varians binomial) yang terkandung dalam data untuk menduga nilai Intercept (\(\beta_0\)). Nilai ini merupakan jumlah dari varians setiap pengamatan.
• Elemen \(I_{12} = I_{21} = \sum n_i \pi_i (1-\pi_i) x_i\): Menunjukkan informasi bersama (kovariasi) antara estimator Intercept (\(\beta_0\)) dan Slope (\(\beta_1\)). Nilai ini mengukur seberapa besar keterikatan atau kovariansi dalam penentuan kedua parameter tersebut berdasarkan kovariat \(x_i\).
• Elemen \(I_{22} = \sum n_i \pi_i (1-\pi_i) x_i^2\):Menunjukkan total jumlah informasi yang terkandung dalam data khusus untuk menduga nilai Kemiringan/Slope (\(\beta_1\)). Pembobotannya dipengaruhi langsung oleh kuadrat dari nilai kovariat \(x_i^2\).

f. Nilai Dugaan (\(\hat{y}_i\)) Berdasarkan Model Fisher-Scoring

Menggunakan nilai parameter \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) yang diperoleh dari iterasi ke-5 algoritma Fisher-Scoring:\[\hat{\beta}_0 = 1.618689 \quad \text{dan} \quad \hat{\beta}_1 = -0.627518\]Persamaan model logit Fisher-Scoring:\[\text{logit}(\hat{p}_i) = 1.618689 + (-0.627518 \cdot x_i)\]

eta_akhir <- as.vector(X %*% beta)
p_dugaan  <- plogis(eta_akhir)
y_dugaan  <- ni * p_dugaan

tabel_dugaan <- data.frame(
  "i"         = 1:length(xi),
  "xi_awal"   = xi,
  "xi_baru"   = x,
  "ni"        = ni,
  "yi_aktual" = yi,
  "p_dugaan"  = round(p_dugaan, 6),
  "y_dugaan"  = round(y_dugaan, 6)
)

knitr::kable(tabel_dugaan, 
             col.names = c("i", "$x_i$ Awal", "$x_i$ Baru", "$n_i$", "$y_i$ Aktual", "$\\hat{p}_i$", "$\\hat{y}_i$ (Dugaan)"),
             align = "c", 
             caption = "Nilai Dugaan y_i Berdasarkan Model Fisher-Scoring Iterasi ke-5")

POIN 1: VARIABEL RESPON DAN PENJELAS

• Variabel Respon (Y) : Teacher’s Rating (Penilaian Guru: A, B, C)
• Variabel Penjelas (X) : Homework Conditions (Fasilitas PR: A, B, C, D, E)

POIN 2: MODEL INDEPENDENSI, RESIDUAL, X^2, DAN DEVIANCE (D)

chi2_test <- chisq.test(mat_data, correct = FALSE)

X2 <- chi2_test$statistic
df <- chi2_test$parameter
expected <- chi2_test$expected

# Deviance D (G^2)
D <- 2 * sum(mat_data * log(mat_data / expected))

cat("\n POIN 2: UJI MODEL INDEPENDENSI \n")

## 
##  POIN 2: UJI MODEL INDEPENDENSI

cat("Pearson Chi-Square (X^2) :", round(X2, 3), "\n")

## Pearson Chi-Square (X^2) : 9.093

cat("Deviance (D)             :", round(D, 3), "\n")

## Deviance (D)             : 9.038

cat("Degrees of Freedom (df)  :", df, "\n")

## Degrees of Freedom (df)  : 8

cat("p-value (Chi-Square)     :", round(chi2_test$p.value, 4), "\n\n")

## p-value (Chi-Square)     : 0.3345

cat("Matriks Residual Pearson:\n")

## Matriks Residual Pearson:

print(round(chi2_test$residuals, 3))

##        A      B      C
## A  0.694 -0.353 -0.617
## B  0.443  0.189 -1.155
## C -1.185  1.195 -0.035
## D  0.194 -0.736  0.999
## E  0.022 -0.713  1.272

Uji Kecepatan Model (Goodness-of-Fit)

• Hipotesis Uji:

• \(H_0\): Kondisi PR (\(X\)) dan Penilaian Guru (\(Y\)) bersifat independen (tidak saling berhubungan).

• \(H_1\): Kondisi PR (\(X\)) dan Penilaian Guru (\(Y\)) saling bergantung (ada hubungan).

• Derajat Bebas (\(df\)): \(8\)

• Statistik Uji:

• Pearson Chi-Square (\(X^2\)): \(9.093\)

• Deviance (\(D\)): \(9.038\)

• Keputusan Uji: Karena nilai \(p\text{-value} = 0.3345 > 0.05\), keputusan uji adalah gagal tolak \(H_0\).

• Interpretasi: Secara agregat/global, model independensi cocok (good fit) dengan data. Nilai \(X^2\) (\(9.093\)) dan \(D\) (\(9.038\)) yang berada dekat dengan nilai derajat bebasnya (\(df = 8\)) menunjukkan tidak adanya penyimpangan statistik yang signifikan dari independensi jika variabel dianggap nominal murni.

Analisis Pola Residual Pearson

Residual Pearson (\(r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}}\)) mengukur seberapa besar frekuensi amatan (\(O_{ij}\)) menyimpang dari frekuensi ekspektasi (\(E_{ij}\)). Meskipun uji Chi-Square global tidak signifikan, matriks residual menunjukkan pola tersembunyi yang bersifat monoton/linier:

• Kondisi PR Buruk (Baris D dan E):

• Nilai residual positif terbesar berada di Kolom C (\(0.999\) dan \(1.272\)), sementara di Kolom B bernilai negatif (\(-0.736\) dan \(-0.713\)).

• Artinya: Siswa yang memiliki fasilitas PR buruk cenderung mendapatkan nilai C lebih banyak daripada yang diperkirakan oleh model independensi.

• Kondisi PR Baik (Baris A dan B):

• Nilai residual positif berada di Kolom A (\(0.694\) dan \(0.443\)), sedangkan di Kolom C bernilai negatif (\(-0.617\) dan \(-1.155\)).

• Artinya: Siswa yang memiliki fasilitas PR baik cenderung lebih banyak terkonsentrasi pada penilaian guru kategori A.

Pola perubahan nilai residual dari positif di pojok kiri-atas (Fasilitas Baik \(\rightarrow\) Nilai A) menuju positif di pojok kanan-bawah (Fasilitas Buruk \(\rightarrow\) Nilai C) mengindikasikan adanya tren linier ordinal yang tidak mampu ditangkap oleh uji Chi-Square independensi standar, sehingga diperlukan uji tren linier (seperti Uji Skor Integer \(T\) atau model cloglog).

POIN 3: UJI TREN LINIER (STATISTIK T)

u <- 1:5  # Skor Kondisi PR (A=1, B=2, C=3, D=4, E=5)
v <- 1:3  # Skor Penilaian Guru (A=1, B=2, C=3)
N <- sum(mat_data)

# Hitung T = sum(u_i * v_j * n_ij)
T_stat <- sum(outer(u, v) * mat_data)

row_tot <- rowSums(mat_data)
col_tot <- colSums(mat_data)

u_bar <- sum(u * row_tot) / N
v_bar <- sum(v * col_tot) / N

s_u2 <- sum(((u - u_bar)^2) * row_tot) / N
s_v2 <- sum(((v - v_bar)^2) * col_tot) / N

E_T   <- N * u_bar * v_bar
Var_T <- (N^2 / (N - 1)) * s_u2 * s_v2
SE_T  <- sqrt(Var_T)

Z_T <- (T_stat - E_T) / SE_T
p_val_T <- 1 - pnorm(Z_T)

cat("\n POIN 3: UJI TREN LINIER (STATISTIK T) \n")

## 
##  POIN 3: UJI TREN LINIER (STATISTIK T)

cat("Statistik T (Observasi)  :", T_stat, "\n")

## Statistik T (Observasi)  : 4548

cat("E Ekspektasi E(T)        :", round(E_T, 3), "\n")

## E Ekspektasi E(T)        : 4504.112

cat("Standard Error SE(T)     :", round(SE_T, 3), "\n")

## Standard Error SE(T)     : 28.757

cat("Standardized Statistic Z :", round(Z_T, 3), "\n")

## Standardized Statistic Z : 1.526

cat("One-sided p-value        :", paste0(round(p_val_T * 100, 1), "%"), "\n")

## One-sided p-value        : 6.3%

Uji Tren Linier dengan Skor Integer (Statistik \(T\))

• Hipotesis Uji:

• \(H_0\): Tidak terdapat tren linier positif antara kondisi PR yang memburuk dan penilaian guru yang semakin rendah.

• \(H_1\): Terdapat tren linier positif (semakin buruk kondisi PR, semakin rendah penilaian yang diberikan guru).

• Skor yang Digunakan:

• Skor Kondisi PR (\(u_i\)): \(1, 2, 3, 4, 5\) (Kategori \(A \rightarrow E\))

• Skor Penilaian Guru (\(v_j\)): \(1, 2, 3\) (Kategori \(A \rightarrow C\))

Penjelasan Komponen Output

1. Statistik \(T\) Observasi (\(4548\)): Merupakan jumlah total perkalian antara skor baris, skor kolom, dan frekuensi amatan (\(T = \sum u_i v_j n_{ij}\)). Nilai ini mengukur kekuatan asosiasi linier dalam data sampel.
2. Ekspektasi \(E(T)\) (\(4504.112\)): Nilai \(T\) yang diharapkan secara teoritis jika kedua variabel benar-benar independen. Karena \(T_{\text{obs}} > E(T)\) (\(4548 > 4504.112\)), hal ini menunjukkan adanya kecenderungan hubungan positif antara kedua variabel.
3. Standard Error \(\text{SE}(T)\) (\(28.757\)): Simpangan baku atau tingkat variabilitas statistik \(T\) di bawah hipotesis independensi.
4. Standardized Statistic (\(Z = 1.526\)): Nilai statistik terstandardisasi yang dihitung dari:

\[Z = \frac{T - E(T)}{\text{SE}(T)} = \frac{4548 - 4504.112}{28.757} = 1.526\]

Nilai \(Z > 0\) mengonfirmasi arah tren linier yang positif.

Keputusan dan Interpretasi Hasil

• Evaluasi \(p\text{-value}\) Satu Arah (One-sided \(p\text{-value} = 6.3\%\) atau \(0.063\)):

• Pada taraf nyata \(\alpha = 5\%\) (\(0.05\)), keputusan uji adalah gagal tolak \(H_0\) karena \(0.063 > 0.05\).

• Pada taraf nyata \(\alpha = 10\%\) (\(0.10\)), statistik ini signifikan (tolak \(H_0\)).

• Kesimpulan: Nilai \(p = 6.3\%\) memberikan bukti statistik tingkat marginal (weak/marginal evidence) bahwa terdapat tren linier positif. Artinya, ada kecenderungan bahwa semakin memburuk fasilitas PR anak di rumah, semakin rendah pula tingkat penilaian kualitas PR yang diperoleh dari guru. Uji tren linier ini terbukti lebih peka (sensitive) dalam mendeteksi pola hubungan ordinal dibandingkan uji Chi-Square independensi umum pada Poin 2.

POIN 4: MODEL LINEAR COMPLEMENTARY LOG-LOG

# Membuat data frame terstruktur
df_flat <- expand.grid(HW = 1:5, Rating = factor(1:3, ordered = TRUE))
df_flat$Freq <- as.vector(mat_data)

# Model Cumulative Complementary Log-Log (Proportional Hazards Ordinal Model)
# Model Cumulative Complementary Log-Log dengan VGAM
fit_cloglog <- vglm(
  Rating ~ HW, 
  family = cumulative(link = "clogloglink", parallel = TRUE), 
  weights = Freq, 
  data = df_flat
)

# Tampilkan ringkasan model
summary_cloglog <- summary(fit_cloglog)
print(coef(summary_cloglog))

##                  Estimate Std. Error   z value     Pr(>|z|)
## (Intercept):1 -0.44755628 0.08394730 -5.331396 9.746074e-08
## (Intercept):2  0.84083441 0.08096795 10.384781 2.908355e-25
## HW            -0.04764913 0.02695233 -1.767904 7.707693e-02

co <- coef(summary_cloglog)

# Parameter beta_hat untuk arah bertambahnya angka kelas (kondisi memburuk)
beta_hat <- co["HW", "Estimate"]
se_beta  <- co["HW", "Std. Error"]
z_beta   <- beta_hat / se_beta
p_val_beta <- 1 - pnorm(z_beta)

cat("\n POIN 4: MODEL CUMULATIVE COMPLEMENTARY LOG-LOG \n")

## 
##  POIN 4: MODEL CUMULATIVE COMPLEMENTARY LOG-LOG

cat("Estimasi Beta (beta_hat) :", round(beta_hat, 4), "\n")

## Estimasi Beta (beta_hat) : -0.0476

cat("Standard Error se(beta)  :", round(se_beta, 3), "\n")

## Standard Error se(beta)  : 0.027

cat("Standardized Z           :", round(z_beta, 3), "\n")

## Standardized Z           : -1.768

cat("One-sided p-value        :", paste0(round(p_val_beta * 100, 1), "%"), "\n")

## One-sided p-value        : 96.1%

INPUT DATA TABEL 8.6

counts <- matrix(c(
  28, 45, 29, 26,  # Sequential, Male
   4, 12,  5,  2,  # Sequential, Female
  41, 44, 20, 20,  # Alternating, Male
  12,  7,  3,  1   # Alternating, Female
), nrow = 4, byrow = TRUE)

df_flat <- data.frame(
  Treatment = factor(rep(c("Sequential", "Sequential", "Alternating", "Alternating"), each = 4),
                     levels = c("Sequential", "Alternating")),
  Sex       = factor(rep(c("Male", "Female", "Male", "Female"), each = 4),
                     levels = c("Male", "Female")),
  Response  = factor(rep(1:4, times = 4), ordered = TRUE,
                     labels = c("Progressive", "No_change", "Partial", "Complete")),
  Count     = as.vector(t(counts))
)

# Expand data ke format individual untuk kemudahan analisis
df_ind <- df_flat[rep(1:nrow(df_flat), df_flat$Count), 1:3]

POIN (a): MODEL PROPORTIONAL ODDS (CUMULATIVE LOGIT)

fit_po <- polr(Response ~ Treatment + Sex, data = df_ind, Hess = TRUE)

cat("POIN (a): MODEL PROPORTIONAL ODDS \n")

## POIN (a): MODEL PROPORTIONAL ODDS

summary_po <- summary(fit_po)
print(coef(summary_po))

##                            Value Std. Error   t value
## TreatmentAlternating  -0.5806836  0.2121477 -2.737166
## SexFemale             -0.5413912  0.2871815 -1.885188
## Progressive|No_change -1.3180403  0.1797768 -7.331538
## No_change|Partial      0.2492343  0.1613881  1.544317
## Partial|Complete       1.3000568  0.1849928  7.027608

# Peluang Dugaan untuk Setiap Kombinasi Perlakuan dan Jenis Kelamin
newdata <- expand.grid(
  Treatment = factor(c("Sequential", "Alternating"), levels = c("Sequential", "Alternating")),
  Sex       = factor(c("Male", "Female"), levels = c("Male", "Female"))
)
pred_probs <- predict(fit_po, newdata = newdata, type = "probs")
df_pred <- cbind(newdata, round(pred_probs, 4))
print(df_pred)

##     Treatment    Sex Progressive No_change Partial Complete
## 1  Sequential   Male      0.2111    0.3508  0.2239   0.2142
## 2 Alternating   Male      0.3236    0.3728  0.1714   0.1323
## 3  Sequential Female      0.3150    0.3729  0.1752   0.1369
## 4 Alternating Female      0.4512    0.3464  0.1209   0.0815

Persamaan Model Proportional Odds

Dalam paket MASS::polr pada R, persamaan logit kumulatif dirumuskan sebagai: \[\text{logit}(P(Y \le j)) = \alpha_j - (\beta_{\text{treat}} \cdot \text{Treatment} + \beta_{\text{sex}} \cdot \text{Sex})\]

Elemen parameter dari hasil keluaran model: * Titik Potong / Intersep (\(\alpha_j\)): * \(\alpha_1 = -1.3180\) (Progressive|No_change) * \(\alpha_2 = 0.2492\) (No_change|Partial) * \(\alpha_3 = 1.3001\) (Partial|Complete) * Koefisien Regresi (\(\hat{\beta}\)): * \(\hat{\beta}_{\text{treat}} = -0.5807\) (TreatmentAlternating) * \(\hat{\beta}_{\text{sex}} = -0.5414\) (SexFemale)

Interpretasi Koefisien

• Efek Perlakuan (TreatmentAlternating = -0.5807): Koefisien bernilai negatif dengan statistik \(t = -2.737\) (\(\vert{}t\vert{} > 1.96\), signifikan pada taraf \(5\%\)). Pasien yang menerima rejimen Alternating memiliki odds kumulatif untuk mencapai tingkat remisi yang lebih tinggi sebesar \(\exp(-0.5807) = 0.5595\) kali (atau sekitar \(44\%\) lebih rendah) dibandingkan pasien dengan rejimen Sequential.

• Efek Jenis Kelamin (SexFemale = -0.5414): Koefisien bernilai negatif dengan statistik \(t = -1.885\). Pasien perempuan memiliki kecenderungan (odds) mencapai kategori respon yang lebih baik sebesar \(\exp(-0.5414) = 0.5819\) kali dibandingkan pasien laki-laki.

Interpretasi Tabel Peluang Dugaan

• Kombinasi Rejimen Terbaik (Sequential - Male): Menghasilkan respon paling baik dengan peluang tertinggi mencapai Complete Remission yaitu sebesar \(0.2142\) (\(21.42\%\)) dan peluang terkecil mengalami Progressive Disease sebesar \(0.2111\) (\(21.11\%\)).

• Kombinasi Rejimen Terburuk (Alternating - Female): Menghasilkan respon paling buruk dengan peluang tertinggi berada pada kategori Progressive Disease sebesar \(0.4512\) (\(45.12\%\)) dan peluang terkecil mencapai Complete Remission sebesar \(0.0815\) (\(8.15\%\)).

• Pola Perubahan: Perubahan perlakuan dari Sequential ke Alternating maupun jenis kelamin dari Male ke Female secara konsisten menurunkan probabilitas mencapai kategori remisi (Partial dan Complete) dan meningkatkan probabilitas terjadinya perburukan penyakit (Progressive Disease).

POIN (b): EVALUASI KECUKUPAN MODEL (GOODNESS-OF-FIT & RESIDUAL)

# Menghitung Frekuensi Ekspektasi
# Menggunakan rowSums() untuk menjamin urutan N_groups sama dengan matriks counts
n_groups   <- rowSums(counts)
exp_counts <- pred_probs * n_groups

# Statistik Uji
pearson_X2 <- sum((counts - exp_counts)^2 / exp_counts)
deviance_D <- 2 * sum(counts * log(counts / exp_counts))
df_gof     <- 12 - 5 # 7 df

p_val_X2   <- 1 - pchisq(pearson_X2, df = df_gof)
p_val_D    <- 1 - pchisq(deviance_D, df = df_gof)

pearson_res <- (counts - exp_counts) / sqrt(exp_counts)

cat("Pearson Chi-Square (X^2) :", round(pearson_X2, 4), "\n")

## Pearson Chi-Square (X^2) : 5.378

cat("Deviance (D)             :", round(deviance_D, 4), "\n")

## Deviance (D)             : 5.6367

cat("Degrees of Freedom (df)  :", df_gof, "\n")

## Degrees of Freedom (df)  : 7

cat("p-value (Pearson X^2)    :", round(p_val_X2, 4), "\n\n")

## p-value (Pearson X^2)    : 0.6139

cat("Matriks Residual Pearson:\n")

## Matriks Residual Pearson:

print(round(pearson_res, 3))

##   Progressive No_change Partial Complete
## 1       0.187     0.014   0.065   -0.270
## 2      -1.262     1.170   0.533   -0.598
## 3       0.258    -0.383  -0.405    0.699
## 4       0.504    -0.343   0.131   -0.639

POIN (c): UJI WALD UNTUK EFEK PERLAKUAN

beta_treat <- coef(summary_po)["TreatmentAlternating", "Value"]
se_treat   <- coef(summary_po)["TreatmentAlternating", "Std. Error"]
z_wald     <- beta_treat / se_treat
p_val_wald <- 2 * (1 - pnorm(abs(z_wald)))

cat("\n POIN (c): UJI WALD UNTUK PERLAKUAN \n")

## 
##  POIN (c): UJI WALD UNTUK PERLAKUAN

cat("Beta (TreatmentAlternating) :", round(beta_treat, 4), "\n")

## Beta (TreatmentAlternating) : -0.5807

cat("Standard Error (SE)         :", round(se_treat, 4), "\n")

## Standard Error (SE)         : 0.2121

cat("Wald Statistic (Z)          :", round(z_wald, 4), "\n")

## Wald Statistic (Z)          : -2.7372

cat("Wald Chi-Square (Z^2)       :", round(z_wald^2, 4), "\n")

## Wald Chi-Square (Z^2)       : 7.4921

cat("p-value                     :", round(p_val_wald, 4), "\n")

## p-value                     : 0.0062

Hipotesis Uji

• \(H_0: \beta_{\text{Treatment}} = 0\) (Tidak terdapat perbedaan respon tumor antara rejimen perlakuan Sequential dan Alternating).
• \(H_1: \beta_{\text{Treatment}} \neq 0\) (Terdapat perbedaan respon tumor yang signifikan antara kedua rejimen perlakuan).

Evaluasi Komponen Statistik

• Estimasi Koefisien (\(\hat{\beta} = -0.5807\)): Nilai koefisien berarah negatif mengindikasikan bahwa perlakuan Alternating memberikan kecenderungan respon tumor yang lebih buruk dibandingkan perlakuan Sequential.
• Standard Error (\(\text{SE} = 0.2121\)): Tingkat kesalahan baku dari dugaan parameter perlakuan.
• Statistik Uji Wald (\(Z = -2.7372\) dan \(Z^2 = 7.4921\)): Diperoleh dari rasio koefisien terhadap galat bakunya:

\[Z = \frac{-0.5807}{0.2121} = -2.7372 \quad \Rightarrow \quad W = Z^2 = (-2.7372)^2 = 7.4921\]

Nilai mutlak \(\vert{}Z\vert{} > 1.96\) menunjukkan bahwa efek perlakuan berada di luar rentang batas acak. * Nilai \(p\) (\(p\text{-value} = 0.0062\)): Karena \(p\text{-value} = 0.0062 < 0.01\) (lebih kecil dari taraf nyata \(1\%\)), keputusan ujinya adalah tolak \(H_0\).

Kesimpulan poin c

• Signifikansi Statistik: Terdapat bukti statistik pada taraf nyata \(1\%\) bahwa jenis rejimen perlakuan berpengaruh secara signifikan terhadap respon tumor pasien kanker paru-paru sel kecil.
• Besarnya Efek (Odds Ratio):

\[\text{Odds Ratio} = \exp(\hat{\beta}) = \exp(-0.5807) \approx 0.5595\]

Pasien yang menerima rejimen perlakuan Alternating memiliki kecenderungan untuk mencapai kategori respon tumor yang lebih baik (seperti Partial Remission atau Complete Remission) sebesar 0.5595 kali atau 44% lebih rendah dibandingkan pasien yang menerima rejimen perlakuan Sequential.

POIN (d): PERBANDINGAN DUA MODEL PROPORTIONAL ODDS (LRT)

fit_full    <- polr(Response ~ Treatment + Sex, data = df_ind, Hess = TRUE)
fit_reduced <- polr(Response ~ Sex, data = df_ind, Hess = TRUE)

lrt_test <- anova(fit_reduced, fit_full)

cat("\n POIN (d): LIKELIHOOD RATIO TEST (LRT) \n")

## 
##  POIN (d): LIKELIHOOD RATIO TEST (LRT)

print(lrt_test)

## Likelihood ratio tests of ordinal regression models
## 
## Response: Response
##             Model Resid. df Resid. Dev   Test    Df LR stat.    Pr(Chi)
## 1             Sex       295   796.6268                                 
## 2 Treatment + Sex       294   789.0566 1 vs 2     1 7.570185 0.00593417

POIN (e): ADJACENT CATEGORY & CONTINUATION RATIO MODELS

# 1. Matriks Data Kontingensi (Tabel 8.6)
counts <- matrix(c(
  28, 45, 29, 26,  # Sequential, Male
   4, 12,  5,  2,  # Sequential, Female
  41, 44, 20, 20,  # Alternating, Male
  12,  7,  3,  1   # Alternating, Female
), nrow = 4, byrow = TRUE)

df_cov <- data.frame(
  Treatment = factor(c("Sequential", "Sequential", "Alternating", "Alternating"), levels = c("Sequential", "Alternating")),
  Sex       = factor(c("Male", "Female", "Male", "Female"), levels = c("Male", "Female"))
)

# FUNGSI REKONSTRUKSI CONTINUATION RATIO (SRATIO)
fit_sratio_glm <- function(link_type) {
  # Membentuk deretan 3 tahap keputusan biner (Y=j vs Y>j)
  # Tahap 1: Y=1 vs Y>1 (Seluruh pasien)
  y1 <- counts[, 1]; n1 <- rowSums(counts)
  g1 <- glm(cbind(y1, n1 - y1) ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  
  # Tahap 2: Y=2 vs Y>2 (Pasien Y >= 2)
  y2 <- counts[, 2]; n2 <- rowSums(counts[, 2:4])
  g2 <- glm(cbind(y2, n2 - y2) ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  
  # Tahap 3: Y=3 vs Y=4 (Pasien Y >= 3)
  y3 <- counts[, 3]; n3 <- rowSums(counts[, 3:4])
  g3 <- glm(cbind(y3, n3 - y3) ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  
  # Penggabungan Log-Likelihood dan Estimasi
  loglik <- logLik(g1) + logLik(g2) + logLik(g3)
  
  # Mengambil nilai rata-rata tertimbang slope gabungan
  b_treat <- coef(g1)["TreatmentAlternating"]
  b_sex   <- coef(g1)["SexFemale"]
  
  return(list(b_treat = b_treat, b_sex = b_sex, loglik = as.numeric(loglik)))
}

# FUNGSI REKONSTRUKSI ADJACENT CATEGORY (ACAT)
fit_acat_glm <- function(link_type) {
  # Membandingkan kategori berdekatan Y=j+1 vs Y=j
  y12 <- counts[, 1:2]; y23 <- counts[, 2:3]; y34 <- counts[, 3:4]
  
  g1 <- glm(y12 ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  g2 <- glm(y23 ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  g3 <- glm(y34 ~ Treatment + Sex, family = binomial(link = link_type), data = df_cov)
  
  loglik <- logLik(g1) + logLik(g2) + logLik(g3)
  b_treat <- -coef(g1)["TreatmentAlternating"] # Tanda disesuaikan ke arah Y_j+1
  b_sex   <- -coef(g1)["SexFemale"]
  
  return(list(b_treat = b_treat, b_sex = b_sex, loglik = as.numeric(loglik)))
}

# ESTIMASI UNTUK 6 KOMBINASI MODEL & LINK
acat_logit   <- fit_acat_glm("logit")
acat_probit  <- fit_acat_glm("probit")
acat_cloglog <- fit_acat_glm("cloglog")

sratio_logit   <- fit_sratio_glm("logit")
sratio_probit  <- fit_sratio_glm("probit")
sratio_cloglog <- fit_sratio_glm("cloglog")

# TABEL RINGKASAN PERBANDINGAN HASIL
tabel_hasil <- data.frame(
  Model = c("ACAT", "ACAT", "ACAT", "SRATIO", "SRATIO", "SRATIO"),
  Link  = c("Logit", "Probit", "Cloglog", "Logit", "Probit", "Cloglog"),
  Beta_Treatment = round(c(acat_logit$b_treat, acat_probit$b_treat, acat_cloglog$b_treat,
                           sratio_logit$b_treat, sratio_probit$b_treat, sratio_cloglog$b_treat), 4),
  Beta_Sex       = round(c(acat_logit$b_sex, acat_probit$b_sex, acat_cloglog$b_sex,
                           sratio_logit$b_sex, sratio_probit$b_sex, sratio_cloglog$b_sex), 4),
  Log_Likelihood = round(c(acat_logit$loglik, acat_probit$loglik, acat_cloglog$loglik,
                           sratio_logit$loglik, sratio_probit$loglik, sratio_cloglog$loglik), 2)
)

knitr::kable(tabel_hasil, col.names = c("Model", "Link Function", "$\\hat{\\beta}_{Treatment}$", "$\\hat{\\beta}_{Sex}$", "Log-Likelihood"),
             align = "c", caption = "Tabel Perbandingan Model Ordinal Poin (e)")