test_that("rf_arvind model fitting, prediction, and diagnostics work for EM and MCMC methods", {
  set.seed(123)
  df <- data.frame(
    y = 10 + rnorm(40),
    x1 = rnorm(40),
    x2 = runif(40)
  )
  
  # EM method
  fit_em <- rf_arvind(
    y ~ x1 + x2,
    data = df,
    method = "em",
    num.trees = 30,
    n_boot = 20,
    seed = 123
  )
  
  expect_s3_class(fit_em, "rf_arvind")
  expect_true(fit_em$theta > 0)
  expect_true(fit_em$mse >= 0)
  expect_true(fit_em$rmse >= 0)
  expect_true(fit_em$mae >= 0)
  expect_true(fit_em$risk >= 0)
  expect_true(is.matrix(fit_em$intervals$theta_ci))
  
  # Prediction with EM fit
  preds_em <- predict(fit_em, newdata = df[1:5, ])
  expect_equal(nrow(preds_em), 5)
  expect_true("PI_90_Lower" %in% colnames(preds_em))
  expect_true("PI_95_Lower" %in% colnames(preds_em))
  expect_true("PI_99_Lower" %in% colnames(preds_em))
  
  # MCMC method
  fit_mcmc <- rf_arvind(
    y ~ x1 + x2,
    data = df,
    method = "mcmc",
    num.trees = 30,
    n_iter = 1000,
    burn_in = 200,
    seed = 123
  )
  
  expect_s3_class(fit_mcmc, "rf_arvind")
  expect_true(fit_mcmc$theta > 0)
  expect_length(fit_mcmc$intervals$hpd_90, 2)
  expect_length(fit_mcmc$intervals$hpd_95, 2)
  expect_length(fit_mcmc$intervals$hpd_99, 2)
  
  # Summary and diagnostics
  sum_em <- summary(fit_em)
  expect_s3_class(sum_em, "summary.rf_arvind")
  expect_true(!is.na(sum_em$ks_stat))
  expect_true(!is.na(sum_em$ad_stat))
  expect_true(!is.na(sum_em$aic))
  expect_true(!is.na(sum_em$bic))
  
  # NA handling check
  df_na <- df
  df_na$y[1] <- NA
  df_na$x1[2] <- NA
  fit_na <- rf_arvind(y ~ x1 + x2, data = df_na, num.trees = 30, n_boot = 10, seed = 123)
  expect_s3_class(fit_na, "rf_arvind")
})