Why does a fixed-effect OLS need unique time elements? The plm function of the plm library in R is giving me grief over having duplicate time-id couples, even when I'm running a model that I don't think should need a time variable at all (see reproducible example below).
I can think of three possibilities:


*

*My understanding of fixed effects regression is wrong, and they really do require unique time indices (or time indices at all!).

*plm() is just being overly-finicky here and should relax this requirement.

*The particular estimation technique that plm() uses--the within transformation--requires time indices, even though the order doesn't seem to matter and the less computationally-efficient version (including dummies in a straight-up OLS model) doesn't need them.


Any thoughts?
set.seed(1)
n <- 1000
test <- data.frame( grp = as.factor(rep( letters, (n/length(letters))+1 ))[seq(n)], x = runif(n), z = runif(n) )
test$y <- with( test, 2*x + 3*z + rnorm(n) )
lm( y ~ x + z, data = test )
lm( y ~ x + z + grp, data = test )

require(plm)
# Model fails if I don't specify a time index, despite effect = "individual"
plm( y ~ x + z, data = test, model = "within", effect="individual", index = "grp" ) 
# Create time variable and add it to the index but still specify individual FE not time FE also
library(plyr)
test <- ddply( test, .(grp), function(dat) transform( dat, t = seq(nrow(dat)) ) )
# Now plm() works; note coefficients clearly include the fixed effects, as they match the lm() version above
plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") ) 
# Scramble time variables and show they don't matter as long as they're unique within a cluster
test <- ddply( test, .(grp), function(dat) transform( dat, t = sample(t) ) )
plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") ) 
# Add a duplicate time entry and show that it causes plm() to fail
test[ 2, "t" ] <- test[ 1, "t" ] 
plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") ) 

Why this matters
I'm trying to bootstrap my model, and when I do the requirement that the index-time pairs be unique is causing headaches which seem unnecessary if (2) is true.
 A: Your understanding of fixed effects regression seems perfectly fine. When you do the within transformation to obtain the fixed effects estimator
$$y_{it} - \overline{y}_{i} = (X_{it} - \overline{X}_i)\beta + \epsilon_{it} - \overline{\epsilon}_i$$
the time-sorting order does not matter because $\overline{y}_{i} = \frac{1}{T}\sum^{T}_{t=1}y_{it}$, $\overline{x}_{i} = \frac{1}{T}\sum^{T}_{t=1}x_{it}$, and $\overline{\epsilon}_{i} = \frac{1}{T}\sum^{T}_{t=1}\epsilon_{it}$ sum out the time component no matter the sorting order within each individual (or firm/country/whatever your $i$ subscript may be).
I'm not an R guy but in Stata you would run into the same problem for duplicate time values in the time variable. Again this wouldn't matter for the fixed effects estimation and in fact you do not even need to specify a time variable.
For example,
webuse nlswork
xtset idcode
xtreg ln_wage age hours, fe

will give you the same estimates as
xtset idcode year
xtreg ln_wage age hours, fe

The sorting order of the time values can sometimes be important for inference though. If you were to use the xtserial command after the above fixed effects regression Stata will tell you
xtserial age
time variable not set, use -tsset varname ...

if you haven't used xtset idcode year before. For this purpose it can be problematic if you have 2 observations for an individual in a given year but you do not know if one observation is dated before/after the other (for instance, if a month or quarter variable is missing).
I'm sure this is not the case for you but sometimes people specify the time variable to be annual when in fact they have monthly data. If they wanted to run such a regression they would need to aggregate the data first to the annual level. Otherwise to solve the duplicate time values problem one would generate a new time variable for year-month combinations. The within estimator itself does not need a specified time component itself.
A: Indeed, plm will not allow you to run a FE model, when there is a lower-level unit (i.e. you want household instead of individual, country instead of states etc). And indeed, there's nothing wrong about doing what you want. 
The trick in this case is just to make the time variable unique, crossing it with the sub-level unit: make a time-individual if you do at household level, or state-year if you do at region level. See similar post: https://stackoverflow.com/questions/43510067/fixed-effects-plm-package-r-multiple-observations-per-year-id/43573731
library(plm)
#> Loading required package: Formula
data("Produc", package = "plm")


Produc$year_state <- paste(Produc$year, Produc$state, sep="_")

## will throw warning
Produc_plm <- pdata.frame(Produc, index = c("region", "year"))
#> Warning in pdata.frame(Produc, index = c("region", "year")): duplicate couples (id-time) in resulting pdata.frame
#>  to find out which, use e.g. table(index(your_pdataframe), useNA = "ifany")

## will throw error:
reg_plm_1 <- plm(gsp ~ pcap, data = Produc_plm)
#> Warning: non-unique values when setting 'row.names': '1-1970', '1-1971',
#> '1-1972', '1-1973', '1-1974', '1-1975', '1-1976', '1-1977', '1-1978',
#> '1-1979', '1-1980', '1-1981', '1-1982', '1-1983', '1-1984', '1-1985',

#> Error in `.rowNamesDF<-`(x, value = value): duplicate 'row.names' are not allowed

Use the trick instead:
Produc_plm2 <- pdata.frame(Produc, index = c("region", "year_state"))
reg_plm_2 <- plm(gsp ~ pcap, data = Produc_plm2)

Let's check with package lfe if we got it right:
library(lfe)
#> Loading required package: Matrix
#> 
#> Attaching package: 'lfe'
#> The following object is masked from 'package:plm':
#> 
#>     sargan
library(broom)
reg_lfe_1 <- felm(gsp ~ pcap|region, data = Produc)
all.equal(as.data.frame(tidy(reg_plm_2)), 
          as.data.frame(tidy(reg_lfe_1)))
#> [1] TRUE

