For small numbers of variables, this problem is easier than it looks: for each group of data you can systematically try all $m!$ permutations of the coefficients and retain the permutation that minimizes the sum of residual squares.
To describe the algorithm in detail, let the data consist of $g$ pairs $X^{(j)}, Y^{(j)}$ where each $X^{(j)}$ is an $n_j$ by $m$ matrix whose columns are the independent values and each $Y^{(j)}$ is an $n_j$ column vector of dependent values. The model is
$$E[Y^{(j)}] = X^{(j)}\sigma_j(\beta)$$
where
- $\beta$ is an $m$-vector of coefficients
- $\sigma_j$ are permutations of the coefficients of $\beta$
- $1 \le j \le g$.
Let $\mathfrak{S}_m$ denote the group of all permutations of the positions of an $m$-vector. To estimate the parameters $\beta$ and $\{\sigma_j\}$, define
$$f_j(\beta) = \min_{\sigma \in \mathfrak{S}_m} ||Y^{(j)} - X^{(j)}\sigma(\beta)||^2.$$
Given $\beta$, this amounts to finding a permutation of $\beta$ which minimizes the residuals in group $j$. One way to find this permutation is to compute $||Y^{(j)} - X^{(j)}\sigma(\beta)||^2$ for each of the $m!$ permutations $\sigma$. This brute force approach is reasonable for small $m$.
The estimate of $\beta$ minimizes the objective $F(\beta) = \sum_{j=1}^g f_j(\beta)$. Because $F$ is smooth almost everywhere, continuous, and quadratic in a neighborhood of is minimum, we may apply a gradient-based numerical solver to find the optimal value $
\widehat\beta$. This can be expected to converge quickly. Thus, the total effort is proportional (roughly) to $g m! \bar n$ where $\bar n = (n_1+\cdots+n_g)/g$ is the average group size.
Having obtained an estimate of $\beta$--which involves estimating the $\sigma_j$ when evaluating the $f_j$--we may then proceed as usual to predict $Y$, compute residuals, perform all diagnostics, etc. A convenient way to implement this would be to permute the rows of the data $X^{(j)}$ once and for all according to the inverse of the estimated $\sigma_j$ and then perform a standard regression on the dataset obtained by combining all the $X^{(j)}, Y^{(j)}$ cases.
The usual least-squares (and maximum likelihood) theory should hold approximately provided the groups are large and not too great in number. (When there is a large number of groups, some of the $\sigma_j$ are likely to be incorrectly estimated, leading to a reduction in the sum of squared residuals: in effect, we should expect to overfit the data.)
As a moderately difficult example I synthesized $1979$ cases with $m=6$ variables (each generated independently from a standard Normal distribution) partitioned into $g=10$ groups averaging almost $200 = \bar n$ per group. I set $\beta=(1,-2,3,-4,5,-6)$ and added iid Normally distributed error $\varepsilon$ with standard deviation $\tau=12$. After $50$ evaluations of $F$, requiring $15$ seconds total, convergence was achieved at the (reasonable) estimate $\widehat\beta = (1.10, -2.15, 2.88, -4.50, 4.98, -6.41)$. The estimate of $\tau$ was $12.12$ (to be compared to the actual realized value of $12.15$). Evidently overfitting was not a problem in this example.
Here is the R code used for the example. It is written in a moderately general way to allow flexible experimentation and even application to real data.
#
# Specify the problem.
#
set.seed(17)
m <- 6 # Number of variables (including any constant)
n.groups <- 10 # Number of groups
n.per.group <- rbinom(n.groups, 250, 4/5) # Group sizes
n <- sum(n.per.group) # Number of cases
eps <- rnorm(n, 0, 2*m) # Error terms
suppressWarnings(beta <- 1:m * c(1,-1)) # Coefficients
f <- t(replicate(n.groups, sample.int(m))) # Permutations
#
# Create data (y, x, g): IV, DVs, and groupings.
#
pred <- function(b, x, group, s) {
#
# Apply the parameters (b, s) to (x, group).
#
unlist(sapply(unique(group), function(g) x[group==g, ] %*% b[s[g, ]]))
}
group <- unlist(sapply(1:n.groups, function(i) rep(i, n.per.group[i])))
x <- matrix(rnorm(m*n), n, m)
dimnames(x) <- list(NULL, paste("x", 1:m, sep="."))
y <- pred(beta, x, group, f) + eps
#
# Define the objective function.
#
permutations <- function(x) {
#
# Returns all the permutations of `x`, one per row.
#
if (length(x) <= 1) return(matrix(x, 1, length(x)))
matrix(sapply(1:length(x), function(i) t(cbind(permutations(x[-i]), x[i]))),
ncol=length(x), byrow=TRUE)
}
fit <- function(y, x, b, sigma=permutations(1:length(b))) {
#
# Returns a permutation `sigma` minimizing |y - x.sigma(b)|^2 along
# with the minimal value.
#
delta <- x %*% apply(sigma, 1, function(s) b[s])
costs <- apply(delta, 2, function(z) {w <- z-y; sum(w*w)})
i <- which.min(costs)
list(permutation=sigma[i,], value=costs[i])
}
cost <- function(y, x, b, group=1, sigma=permutations(1:length(b))) {
#
# Returns the total cost and the array of permutations minimizing it
# for IVs `x`, DV `y`, and coefficients `b`. The grouping of `x`
# is given by `group`. The feasible permutations are rows of `sigma`.
#
sapply(unique(group), function(g) {
i <- group==g
fit(y[i], x[i, ], b, sigma)
})
# Returns one column per group.
}
objective <- function(b, y, x, g, sigma=permutations(1:length(b))) {
#
# Returns the value for coefficients `b` and IVs `x`, DV `y`, groups `g`.
#
sum(unlist(cost(y, x, b, g, sigma)["value", ]))
}
#
# Find *some* starting value (albeit a poor one).
#
b <- coef(lm(y ~ x - 1))
b <- beta
#
# Optimize.
#
sigma <- permutations(1:length(b))
system.time(model <- optim(b, function(b) objective(b, y, x, group, sigma),
method="BFGS"))
#
# The fitted coefficients should be compared to `beta` *up to a permutation*!
#
model
sqrt(model$value / (n-1)) # (Under)estimates the SD of `eps` $
sd(eps)
#
# Diagnostics and checks.
#
f.hat <- t(sapply(cost(y, x, model$par, group, sigma)["permutation", ], identity))
y.hat <- pred(model$par, x, group, f.hat)
pairs(cbind(x, y, y.hat), pch=".")
#
# Additional checks.
#
cost(y, x, b, group, sigma)
objective(beta, y, x, group, sigma)
objective(b, y, x, group, sigma)
objective(model$par, y, x, group, sigma)
model$value
