R Under development (unstable) (2024-02-21 r85967 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(lfe)
Loading required package: Matrix
> set.seed(43)
> options(lfe.threads=2,digits=5,warn=1)
> 
> g1 <- 80
> g2 <- 20
> g3 <- 12
> N <- 1000
> clu1 <- sample(g1,N, replace=TRUE)
> clu2 <- (clu1 + sample(7,N,replace=TRUE)-1) %% g2
> clu3 <- (clu2 + sample(3,N,replace=TRUE)-1) %% g3
> clu1 <- factor(clu1)
> clu2 <- factor(clu2)
> clu3 <- factor(clu3)
> # group specific covariate effects
> ceff1 <- rnorm(nlevels(clu1), sd=0.5)[clu1]
> ceff2 <- rnorm(nlevels(clu2), sd=0.4)[clu2]
> ceff3 <- rnorm(nlevels(clu3))[clu3]
> 
>  # group specific errors
> err1 <- rnorm(nlevels(clu1), sd=0.8)[clu1]
> err2 <- rnorm(nlevels(clu2))[clu2]
> err3 <- rnorm(nlevels(clu3), sd=0.5)[clu3]
> 
> x1 <- ceff1 + 0.3*ceff2 + rnorm(N)
> x2 <- ceff2 + 0.2*ceff3 + rnorm(N)
> x3e <- ceff3 + 0.2*(ceff2+ceff1) + rnorm(N)
> 
> f1 <- factor(sample(8,N,replace=TRUE))
> x3 <- as.vector(as(f1,'sparseMatrix') %*% x3e)[f1]/tabulate(f1)[f1]
> err <- err1 + err2 + err3 + abs(x1+x2*x3)*rnorm(N)
> y <- x1 + x2 + x3 + err
> data <- data.frame(y,x1,x2,x3,f1,clu1,clu2,clu3)
> clu <- list('clu1', 'clu2', 'clu3')
> summary(felm(y ~ x1 + x2 + f1|0|0|clu1+clu2+clu3, data))
Warning in newols(mm, nostats = nostats[1], exactDOF = exactDOF, onlyse = onlyse,  :
  Negative eigenvalues set to zero in multiway clustered variance matrix. See felm(...,psdef=FALSE)
Warning in chol.default(mat, pivot = TRUE, tol = tol) :
  the matrix is either rank-deficient or not positive definite

Call:
   felm(formula = y ~ x1 + x2 + f1 | 0 | 0 | clu1 + clu2 + clu3,      data = data) 

Residuals:
   Min     1Q Median     3Q    Max 
-6.408 -1.298  0.066  1.323  8.091 

Coefficients:
            Estimate Cluster s.e. t value Pr(>|t|)    
(Intercept)   0.3142       0.4159    0.76     0.45    
x1            1.0244       0.0789   12.98   <2e-16 ***
x2            0.8472       0.1173    7.22    1e-12 ***
f12           0.3737       0.2062    1.81     0.07 .  
f13           0.2362       0.1796    1.32     0.19    
f14          -0.2908       0.3184   -0.91     0.36    
f15           0.2397       0.2298    1.04     0.30    
f16          -0.3411       0.2236   -1.53     0.13    
f17           0.0382       0.2242    0.17     0.86    
f18           0.0986       0.1340    0.74     0.46    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.97 on 990 degrees of freedom
Multiple R-squared(full model): 0.372   Adjusted R-squared: 0.366 
Multiple R-squared(proj model): 0.372   Adjusted R-squared: 0.366 
F-statistic(full model, *iid*):65.1 on 9 and 990 DF, p-value: <2e-16 
F-statistic(proj model): 1.31e+04 on 9 and 11 DF, p-value: <2e-16 


> #gclu <- structure(clu, method='gaure')
> #summary(felm(y ~ x1 + x2 + f1|0|0|clu1+clu2+clu3, data, cmeth='gaure'))
> #summary(est <- felm(y ~ x1 + x2 | f1, data, clustervar=gclu))
> #ef <- structure(function(x,addnames) {
> #    c(x[1],x[2:8]-x[1])
> #}, verified=TRUE)
> #getfe(est,ef=ef,se=TRUE, bN=200)
> 
> proc.time()
   user  system elapsed 
   1.10    0.09    1.20