diff --git a/R/fitstats.R b/R/fitstats.R
index e6ee23b8..4a977f26 100644
--- a/R/fitstats.R
+++ b/R/fitstats.R
@@ -1619,12 +1619,23 @@ kp_stat = function(x){
   # internal function => x must be a fixest object
   #
   # The code here is a translation of the ranktest.jl function from the Vcov.jl package
-  # from @matthieugomez (see https://github.com/matthieugomez/Vcov.jl)
-  #
+  # from @matthieugomez (see https://github.com/FixedEffects/Vcov.jl)
 
 
   if(!isTRUE(x$iv) || !x$iv_stage == 2) return(NA)
 
+  n_endo = length(x$iv_first_stage)
+  n_inst = x$iv_n_inst
+  VCOV_TYPE = attr(x$cov.scaled, "type")
+  if (!identical(VCOV_TYPE, "IID") && n_inst != n_endo) {
+    warning(paste0(
+      "KP calculation is only implemented for the cases:\n",
+      " * vcov type is IID or\n",
+      " * number of endogenous variables == number of instruments"
+    ))
+    return(NA)
+  }
+
   # Necessary data
 
   X_proj = as.matrix(resid(summary(x, stage = 1)))
@@ -1637,8 +1648,8 @@ kp_stat = function(x){
   Z = model.matrix(x, type = "iv.inst")
   Z_proj = proj_on_U(x, Z)
 
-  k = n_endo = ncol(X_proj)
-  l = n_inst = ncol(Z)
+  k = n_endo
+  l = n_inst
 
   # We assume l >= k
   q = min(k, l) - 1
@@ -1650,7 +1661,7 @@ kp_stat = function(x){
   } else {
     PI = coef(summary(x, stage = 1))
   }
-  PI = PI[, colnames(PI) %in% x$iv_inst_names_xpd, drop = FALSE]
+  PI = as.matrix(PI[, colnames(PI) %in% x$iv_inst_names_xpd, drop = FALSE])
 
   Fmat = chol(crossprod(Z_proj))
   Gmat = chol(crossprod(X_proj))
@@ -1691,8 +1702,6 @@ kp_stat = function(x){
   # There is need to compute the vcov specifically for this case
   # We do it the same way as it was for x
 
-  VCOV_TYPE = attr(x$cov.scaled, "type")
-
   if(identical(VCOV_TYPE, "IID")){
     vlab = chol(tcrossprod(kronv) / nrow(X_proj))
 
diff --git a/tests/fixest_tests.R b/tests/fixest_tests.R
index cd46b457..3077fb76 100644
--- a/tests/fixest_tests.R
+++ b/tests/fixest_tests.R
@@ -2966,11 +2966,16 @@ names(base) = c("y", "x1", "x_endo_1", "x_inst_1", "fe")
 set.seed(2)
 base$x_inst_2 = 0.2 * base$y + 0.2 * base$x_endo_1 + rnorm(150, sd = 0.5)
 base$x_endo_2 = 0.2 * base$y - 0.2 * base$x_inst_1 + rnorm(150, sd = 0.5)
+base$cl_var = rep(1:50, 3) # using `fe` to cluster causes a vcov PSD warning
 
 # Checking a basic estimation
 est_iv = feols(y ~ x1 | x_endo_1 + x_endo_2 ~ x_inst_1 + x_inst_2, base)
 
-fitstat(est_iv, ~ f + ivf + ivf2 + wald + ivwald + ivwald2 + wh + sargan + rmse + g + n + ll + sq.cor + r2)
+fitstat(est_iv, ~ f + ivf + ivf2 + wald + ivwald + ivwald2 + wh + sargan + kpr + rmse + g + n + ll + sq.cor + r2)
+test(c(length(est_iv$iv_first_stage), est_iv$iv_n_inst), c(2, 2)) # kpr relies on accessing these
+
+est_iv_uneven = feols(y ~ x1 | x_endo_1 ~ x_inst_1 + x_inst_2, base, cluster = "cl_var")
+fitstat(est_iv_uneven, ~ f + ivf + ivf2 + wald + ivwald + ivwald2 + wh + sargan + kpr + rmse + g + n + ll + sq.cor + r2)
 
 est_fe = feols(y ~ x1 | fe, base)