correct bug with fixed effects extraction

R/fepoisson.R

@@ -23,3 +23,27 @@ fepoisson <- function(
     beta.start = beta.start, eta.start = eta.start, control = control
   )
 }
+
+# fequasipoisson <- function(
+#     formula = NULL,
+#     data = NULL,
+#     weights = NULL,
+#     beta.start = NULL,
+#     eta.start = NULL,
+#     control = NULL) {
+#   # Fit the model using standard Poisson assumptions
+#   fit <- feglm(
+#     formula = formula, data = data, weights = weights, family = poisson(),
+#     beta.start = beta.start, eta.start = eta.start, control = control
+#   )
+
+#   # Estimate the dispersion parameter (phi)
+#   fitted_values <- predict(object, type = "response")
+#   residuals <- unlist(object$data[, 1], use.names = FALSE) - fitted_values
+#   phi <- sum((residuals^2) / fitted_values) / fit$df.residual?
+
+#   # Adjust model diagnostics for Quasi Poisson
+#   fit$std.errors <- sqrt(phi) * fit$std.errors
+
+#   return(fit)
+# }

dev/fixedeffects.r

@@ -0,0 +1,46 @@
+load_all()
+
+object <- fepoisson(cyl ~ mpg | am, data = mtcars)
+
+alpha.tol <- 1.0e-08
+
+if (is.null(object)) {
+  stop("'object' has to be specified.", call. = FALSE)
+} else if (isFALSE(inherits(object, "felm")) &&
+  isFALSE(inherits(object, "feglm"))) {
+  stop(
+    "'fixed_effects' called on a non-'felm' or non-'feglm' object.",
+    call. = FALSE
+  )
+}
+
+# Extract required quantities from result list
+beta <- object[["coefficients"]]
+data <- object[["data"]]
+formula <- object[["formula"]]
+lvls.k <- object[["lvls.k"]]
+nms.fe <- object[["nms.fe"]]
+k.vars <- names(lvls.k)
+k <- length(lvls.k)
+eta <- object[["eta"]]
+
+# Extract regressor matrix
+X <- model.matrix(formula, data, rhs = 1L)[, -1L, drop = FALSE]
+nms.sp <- attr(X, "dimnames")[[2L]]
+attr(X, "dimnames") <- NULL
+
+# Generate auxiliary list of indexes for different sub panels
+k.list <- get_index_list_(k.vars, data)
+
+# Recover fixed effects by alternating the solutions of normal equations
+pie <- eta - solve_y_(X, beta)
+
+fe.list <- as.list(get_alpha_(pie, k.list, alpha.tol))
+
+for (i in seq.int(k)) {
+  fe.list[[i]] <- as.vector(fe.list[[i]])
+  names(fe.list[[i]]) <- nms.fe[[i]]
+}
+names(fe.list) <- k.vars
+
+fe.list

src/01_center_variables.cpp

@@ -18,7 +18,7 @@ center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
   const int N = V.n_rows;
   const int P = V.n_cols;
   const int K = klist.size();
-  const double sw = arma::accu(w);
+  const double sw = accu(w);
 
   // Auxiliary variables (storage)
   double delta, denom, meanj, num, wt;
@@ -74,7 +74,7 @@ center_variables_(const doubles_matrix<> &V_r, const doubles &v_sum_r,
       }
 
       // Check convergence
-      delta = arma::accu(arma::abs(x - x0) / (1.0 + arma::abs(x0)) % w) / sw;
+      delta = accu(abs(x - x0) / (1.0 + abs(x0)) % w) / sw;
       if (delta < tol) {
         break;
       }

src/02_get_alpha.cpp

@@ -15,17 +15,17 @@
   Mat<double> y(N, 1);
 
   // Generate starting guess
-  arma::field<arma::vec> Alpha(K);
+  field<Col<double>> Alpha(K);
   for (k = 0; k < K; k++) {
     J = as_cpp<list>(klist[k]).size();
-    Alpha(k) = arma::zeros(J);
+    Alpha(k) = zeros(J);
   }
 
   // Start alternating between normal equations
-  arma::field<arma::vec> Alpha0(arma::size(Alpha));
+  field<Col<double>> Alpha0(size(Alpha));
 
   // Create vector for alpha
-  arma::vec alpha(J);
+  Col<double> alpha(J);
 
   int interruptCheckCounter = 0;
 
@@ -81,8 +81,8 @@
     num = 0.0;
     denom = 0.0;
     for (k = 0; k < K; k++) {
-      num += arma::accu(arma::pow(Alpha(k) - Alpha0(k), 2));
-      denom += arma::accu(arma::pow(Alpha0(k), 2));
+      num += accu(pow(Alpha(k) - Alpha0(k), 2));
+      denom += accu(pow(Alpha0(k), 2));
     }
     crit = sqrt(num / denom);
     if (crit < tol) {
@@ -93,7 +93,9 @@
   // Return alpha
   writable::list Alpha_r(K);
   for (k = 0; k < K; k++) {
-    Alpha_r[k] = as_doubles_matrix(Alpha(k));
+    // Alpha_r[k] = as_doubles_matrix(Alpha(k));
+    Alpha_r[k] = as_doubles(Alpha(k));
   }
+
   return Alpha_r;
 }