From 7b4fe3b664176849a2a17d26fd3af122641f0cec Mon Sep 17 00:00:00 2001
From: Yiqing Xu <yiqingxu@ucsd.edu>
Date: Sun, 12 Feb 2017 23:05:15 -0800
Subject: [PATCH] Fixing bugs

When the sum of weights equals 0 in certain sub-groups.
---
 R/interflex.R   | 55 +++++++++++++++++++++++++++----------------------
 src/fastplm.cpp | 45 ++++++++++++++++++++++++----------------
 2 files changed, 57 insertions(+), 43 deletions(-)

diff --git a/R/interflex.R b/R/interflex.R
index 5419018..729db43 100644
--- a/R/interflex.R
+++ b/R/interflex.R
@@ -422,7 +422,8 @@ inter.binning<-function(
                         xlim = NULL,
                         ylim = NULL,
                         interval = NULL,
-                        Xdistr = "histogram" # c("density","histogram")
+                        Xdistr = "histogram", # c("density","histogram")
+                        wald = TRUE
                         ){
 
 
@@ -1061,9 +1062,11 @@ inter.binning<-function(
                 groupT=data.aug[,time],
                 pairwise=pairwise)$vcov
     }
-    library(lmtest)
-    wald <- waldtest(mod.re, mod.un, test="Chisq", vcov=v)
-    p.wald <- round(wald[[4]][2],4)
+    if (wald == TRUE) {
+        library(lmtest)
+        wald <- waldtest(mod.re, mod.un, test="Chisq", vcov=v)
+        p.wald <- round(wald[[4]][2],4)
+    }
      
     ##################################
     ## storage
@@ -1078,7 +1081,9 @@ inter.binning<-function(
     if (nbins%in%c(2,3)) {
         out<-c(out,list(p.twosided=p.twosided)) 
     }
-    out <- c(out, list(p.wald = p.wald))
+    if (wald == TRUE) {
+        out <- c(out, list(p.wald = p.wald))
+    }
     if (figure==TRUE) {out<-c(out,list(graph=p1))} 
     return(out)
 }
@@ -1629,27 +1634,27 @@ coefs <- function(data,bw,Y,X,D,
     result <- matrix(NA, length(X.eval), (5 + length(Z)))
     colnames(result)<-c("X","Intercept","ME","x","Dx",Z)
     
-    ## formula
-    if (noZ) { ## no covariates
-        formula <- as.formula("y ~ d + xx + d_xx")
-    } else { ## with covariates
-        formula<- as.formula(paste("y ~ d + xx + d_xx + ",paste(Z,collapse="+")))               
-    } 
-    ## weighted least squares without fixed effect
-    wls<-function(x, dat, weights){
-        dat1<-dat
-        dat1$xx<-dat1$x-x
-        dat1$d_xx <- dat1$d * dat1$xx
-        dat1$w <- dnorm(dat1$xx/bw) * weights 
-        reg <- lm(formula, data=dat1, weights = w)
-        result <- c(x, reg$coef)
-        result[which(is.na(result))] <- 0 
-        return(result)   
-    } ## end of WLS function       
-
     ## main algorithm
     if (noFE) { # no fixed effects, wls
-         ## data to be used
+        ## formula
+        if (noZ) { ## no covariates
+            formula <- as.formula("y ~ d + xx + d_xx")
+        } else { ## with covariates
+            formula<- as.formula(paste("y ~ d + xx + d_xx + ",paste(Z,collapse="+")))               
+        } 
+        ## weighted least squares without fixed effect
+        wls<-function(x, dat, weights){
+            dat1<-dat
+            dat1$xx<-dat1$x-x
+            dat1$d_xx <- dat1$d * dat1$xx
+            dat1$w <- dnorm(dat1$xx/bw) * weights 
+            reg <- lm(formula, data=dat1, weights = w)
+            result <- c(x, reg$coef)
+            result[which(is.na(result))] <- 0 
+            return(result)   
+        } ## end of WLS function
+        
+        ## data to be used
         dat<-data[, c(Y, D, X, Z)]
         colnames(dat)[1:3] <- c("y","d","x")
         ## estimation (a loop)
@@ -1659,7 +1664,7 @@ coefs <- function(data,bw,Y,X,D,
         result <- data.frame(result)    
     } else{ # with fixed effects
         ## data to be used
-        dat <- data[, c(Y, D, X, Z)]
+        dat <- data[, c(Y, D, X, Z)] 
         dat.FE <- as.matrix(data[,FE],drop = FALSE)
         for (i in 1:length(X.eval)) {
             xx <- dat[,X]-X.eval[i]
diff --git a/src/fastplm.cpp b/src/fastplm.cpp
index 7ab53e4..fe7b19b 100644
--- a/src/fastplm.cpp
+++ b/src/fastplm.cpp
@@ -7,7 +7,7 @@ using namespace Rcpp;
 List fastplm(arma::mat data,
              arma::mat FE,
              arma::colvec weight,
-             int FEcoefs = 0
+             int FEcoefs = 0 
              ){
       
   // parse data
@@ -23,16 +23,16 @@ List fastplm(arma::mat data,
   /*if(weight.is_empty()){
     weight = arma::ones(n, 1);
   }*/
-
+  
   // declare variables
   arma::colvec y;  // n*1
   arma::mat X; // n*p
   arma::colvec resid; // n*1
   arma::colvec e; // n*1 (with fixed effects)
   arma::colvec coeff; // coefficient (full)
-  arma::colvec se; // SE (full)
+  //arma::colvec se; // SE (full)
   arma::colvec coef; // coefficient
-  arma::colvec stderror; // SE
+  //arma::colvec stderror; // SE
   int df; // degrees of freedom
   double sig2; // sigma2
   double mu; // grand mean
@@ -40,7 +40,7 @@ List fastplm(arma::mat data,
   arma::mat W; // big weighting matrix to calculate fixed effects
   arma::colvec alphas; // fixed effect coefficients
   
-  // count total number of groups (loss of degrees of freedom)
+  /* count total number of groups (loss of degrees of freedom) */
   arma::mat FEvalues;  
   for(int ii=0; ii<m; ii++){ // cluster
     arma::colvec fe = arma::unique(FE.col(ii));   
@@ -52,7 +52,8 @@ List fastplm(arma::mat data,
   }
   int gtot = FEvalues.n_rows;
 
-  // FWL-MAP iteration
+  
+  /* FWL-MAP iteration */
   int niter = 0;
   while ((diff > 1e-5) & (niter<50)) { 
 
@@ -81,7 +82,11 @@ List fastplm(arma::mat data,
       // take average
       for(int i=0; i<(p+1); i++){
         for(int j=0; j<g; j++){
-          mean_value(j,i)=mean_value(j,i)/fe_length(j); 
+           if (fe_length(j)!=0) {
+            mean_value(j,i)=mean_value(j,i)/fe_length(j);
+          } else{
+            mean_value(j,i)=0;
+          } 
         }
       }
       // demean
@@ -101,27 +106,29 @@ List fastplm(arma::mat data,
     data_old = data;
     niter++;
   }
-      
-  // Estimation
+
+ 
+ 
+  /* Estimation */
   y = data.col(0);  // n*1
   if (p>0) {
     X = data.cols(1, p); // n*p
     //store coefficents and check variation of X
     coeff =arma::zeros(p,1); //store coefficients
-    se =arma::zeros(p,1);  
+    //se =arma::zeros(p,1);
+
     // check X variation
     int cc = p;
     for(int i=0;i<cc;i++){
       arma::colvec var =arma::unique(X.col(i));
       if(var.n_rows==1){ // no variation
         coeff(i)=arma::datum::nan ;
-        se(i)=arma::datum::nan ;
+        //se(i)=arma::datum::nan ;
         X.shed_col(i);
         i=i-1;
         cc=cc-1;
       }
     }
-    
     // weighted
     y = y%sqrt(weight);
     for(int i=0;i<cc;i++){
@@ -136,21 +143,23 @@ List fastplm(arma::mat data,
   else {
     resid = y;
   }
+ 
       
   // std.err.
   df = n - gtot - p;
   sig2 = arma::as_scalar(resid.t()*resid/df);
-  if (p>0) { 
-    stderror = arma::sqrt(sig2 * arma::diagvec(arma::inv(arma::trans(X)*X))); 
-  } 
+  // if (p>0) { 
+  //   stderror = arma::sqrt(sig2 * arma::diagvec(arma::inv(arma::trans(X)*X))); 
+  // } 
       
   // fill in non-NaN coefficients
   int count=0;
   //String label = xname(0); 
   for(int i=0; i<p; i++){
-    if(coeff(i)==0||se(i)==0){
+    if(coeff(i)==0){
+    //  if(coeff(i)==0||se(i)==0){
       coeff(i)=coef(count);
-      se(i)=stderror(count); 
+      //  se(i)=stderror(count); 
       count=count+1;
     }
   }
@@ -217,7 +226,7 @@ List fastplm(arma::mat data,
   List output;
   if (p > 0) {
     output["coefficients"] = coeff;
-    output["stderr"] = se;
+    //output["stderr"] = se;
   }
   output["residuals"] = resid;
   output["niter"] = niter;