From 8a499ae8bca8dd2108f4115768fd8797039f46df Mon Sep 17 00:00:00 2001
From: Marder <fm58hufi@usr.idiv.de>
Date: Tue, 8 Oct 2024 15:03:48 +0200
Subject: [PATCH] funtion test coef diff
---
Scripts/compare_coef.R | 38 ++++++++++++++++++++++++++++++++++++++
1 file changed, 38 insertions(+)
create mode 100644 Scripts/compare_coef.R
diff --git a/Scripts/compare_coef.R b/Scripts/compare_coef.R
new file mode 100644
index 0000000..cd844db
--- /dev/null
+++ b/Scripts/compare_coef.R
@@ -0,0 +1,38 @@
+
+#load model
+model <- readRDS("Estimation_results/mxl/prediction/MXL_wtp_Prediction matching all cp_model.rds")
+
+compare_coefficients <- function(model, coef_name1, coef_name2) {
+ # Extract coefficients and standard errors
+ coef <- summary(model)$estimate
+ beta1 <- coef[coef_name1, "Estimate"]
+ beta2 <- coef[coef_name2, "Estimate"]
+
+
+ # Assuming model$robse contains robust standard errors
+ se1 <- model$se[coef_name1]
+ se2 <- model$se[coef_name2]
+ # Assuming model$robse contains robust standard errors
+ robse1 <- model$robse[coef_name1]
+ robse2 <- model$robse[coef_name2]
+ #Caluclate Difference
+ difference <- beta1-beta2
+ # Calculate Z-score
+ z_value <- (beta1 - beta2) / sqrt(se1^2 + se2^2)
+
+ # Calculate p-value
+ p_value <- 2 * (1 - pnorm(abs(z_value)))
+ # Calculate Z-score
+ rob_z_value <-(beta1 - beta2) / sqrt(robse1^2 + robse2^2)
+
+ # Calculate p-value
+ rob_p_value <- 2* (1 - pnorm(abs(rob_z_value)))
+
+ # Return results as a list
+ return(list(difference = difference ,z_value = z_value, p_value = p_value ,rob_z_value = rob_z_value, rob_p_value = rob_p_value))
+}
+
+# Example usage
+compare_coefficients(model, "mu_nat_vol_treated", "mu_nat_treat_pred")
+compare_coefficients(model, "mu_nat_no_info", "mu_nat_treat_not_pred")
+
--
GitLab