# Manual two-way demeaning panel data in R

I'm trying to manually demean panel data by both time demeaning and cross-sectional demeaning, yet haven't been able to do it well. The best thing I can think of is looping through the means for each year and then for each entity, plus for each variable. Which seems like it would get slow with a lot of a data (especially if I plan on doing simulations). Is there a better way to do this?

I couldn't find a function in R to do it either.

Here is what I have so far: (DATA is a data frame that of the format (n,T, Variables) where n is the entity id and T is the time id);

demean=function(DATA){

names=colnames(DATA)
T=(max(DATA[,2])-min(DATA[,2])+1)
N=max(DATA[,1])

##Cross-Sectional Demeaning
widedata=reshape(DATA, direction="wide", v.names=names[-c(1:2)],     idvar=names[2], timevar=names[1])
crossmean=matrix(NA, ncol=length(colnames(widedata))-1)
crossmean[,1:length(t(crossmean[1,]))]=colMeans(widedata[,-c(1)], na.rm=TRUE)
crosswidedata=widedata
for(i in 1:T){
crosswidedata[i,-c(1)]=widedata[i,-c(1)]-crossmean
}

crossdemeaned=reshape(crosswidedata, direction="long", times=names[2] )

##Time Demeaning
widedata=reshape(crossdemeaned, direction="wide", v.names=names[-c(1:2)], idvar=names[1], timevar=names[2])
timemean=matrix(NA, ncol=length(colnames(widedata))-1)
timemean[,1:length((t(timemean[1,])))]=colMeans(widedata[,-c(1)], na.rm=TRUE)
timewidedata=widedata
for(i in 1:N){
timewidedata[i,-c(1)]=widedata[i,-c(1)]-timemean
}

demeaned=reshape(timewidedata, direction="long", times=names[2] )

demeaned=demeaned[order(demeaned[,names[2]],demeaned[,names[1]]),]

return(demeaned)
}


I simulate a random spatial autoregressive panel data set with fixed and time effects, demean the data, and then run a test within regression of the demeaned data as following:

    ##Creating a Random Spatial Autoregressive Panel Dataset with Fixed and Time Effects

library(spdep)
library(splm)

set.seed(44222)

################################################################
# RANDOMLY INSERT A CERTAIN PROPORTION OF NAs INTO A DATAFRAME #
################################################################
NAins <-  NAinsert <- function(df, prop = .1){
n <- nrow(df)
m <- ncol(df)
num.to.na <- ceiling(prop*n*m)
id <- sample(0:(m*n-1), num.to.na, replace = FALSE)
rows <- id %/% m + 1
cols <- id %% m + 1
sapply(seq(num.to.na), function(x){
df[rows[x], cols[x]] <<- NA
}
)
return(df)
}
############## df means data frame, and prop is the proportion of missing ##data desired

############################
########CREATES A RANDOM WEIGHT MATRIX
############################
DAG.random <- function(v, nedges=1) {
edges.max <- v*(v-1)/2
# Assert length(v)==1 && 1 <= v
# Assert 0 <= nedges <= edges.max
index.edges <- lapply(list(1:(v-1)), function(k) rep(k*(k+1)/2, v-k))
index.edges <- index.edges[[1]] + 1:edges.max
}

###################################
###### Create Data ################
###################################

n=50                                    ##Number of Observations
T=20                                    ##Number of years
connect=6                               ##average number of connections
W=DAG.random(n,n*.5*connect)            ##Creates a random weight matrix
W=W+t(W)
W=sweep(W, 1, rowSums(W), FUN="/")      ##Row Normalizes the Random Weight

v=3                                     ##Number of Independent Variables

X=matrix(data=NA, nrow=n*T, ncol=(v+2))
for(jjj in 1:T){
X[((jjj-1)*n+1):(jjj*n),1]=c(1:n)
X[((jjj-1)*n+1):(jjj*n),2]=jjj
for(jj in 3:(v+2)){
X[((jjj-1)*n+1):(jjj*n),jj]=rnorm(n, 1, 2)
}
}

sigma=1                             ##standard deviation of model
p=.9                                ##coefficient on spatial lag of Y
B0=matrix(1, nrow=n, ncol=1)
B1=matrix(2, v , 1)
B2=matrix(.5, v, 1)
yrcoef=matrix(runif(1:(T), -20,20), (T), 1)     #Year Fixed Effects
yrcoef[1]=0
alpha=runif(n, -20,20)                      # Nation Fixed Effects

e=matrix(rnorm(n*T,0,sigma), n*T, 1)

inv=solve((diag(n)-p*W))
Y=matrix(data=NA, nrow=n*T, ncol=1)

for(jjj in 1:T){
invY= as.matrix(alpha) + X[((jjj-1)*n+1):(jjj*n),(3:(2+v))] %*% B1 +    yrcoef[jjj,1] + e[((jjj-1)*n+1):(jjj*n),1]
Y[((jjj-1)*n+1):(jjj*n),1]=inv %*% invY
}

Y=data.frame(Y)
X=data.frame(X)
DATA=cbind(X,Y)
names=c("N", "T", "X1", "X2", "X3", "Y")
colnames(DATA)=names

DATA=demean(DATA)

model1=as.formula(paste(names[length(names)], paste(names[-c(1:2,length(names))], collapse= " + "), sep=" ~ "))

lw=mat2listw(W)

mod=spml(model1, data=DATA, listw=lw, model="within", effect=c("twoways"),    lag=TRUE, spatial.error="none")

effects(mod)


However, the time period effects are non-zero (the spatial fixed effects are numerically zero). I can't seem to figure out why the spatial fixed effects is working while the time period effects portion is not.

Any thoughts?

• I'm voting to close this question as off-topic because it is about how to do something in R. (There is a kind of question about whether you can "demean" in some other way, but I don't know what that is aiming at. Means are built-in functions in any decent software, so explicit looping should not be needed.) Without a reproducible example or an attempt at code, this is not especially suited for migration to Stack Overflow. Please review advice on software-related questions in the Help Center. Feb 3 '15 at 16:08
• @Nick This actually looks to me like an inherently statistical problem disguised as a 'how do I do this in language X' question; if this wasn't closed I'd probably have posted something about using regression/ANOVA to remove the means, which I imagine ought to be possible whatever platform one works in. Feb 3 '15 at 23:44
• @Glen_b The question was edited by adding large chunks of code and then asking why the results are puzzling. The statistical question you hope is buried inside may be more obvious to the many active members who read R, but it still looks more like a computing question to me. Feb 5 '15 at 10:02
• @Nick The edits do seem to make it more clearly off topic Feb 5 '15 at 11:09