# R/rem.R In bpm: Basic functions for panel data analysis

```"rem" <-
function(formula,id,data,...)
#Random-effects GLS regression
{
cl <- cl2 <-cl3 <- match.call()
#within effect model
cl2[[1]]<-as.name("fem")
fix.eff <- eval(cl2,parent.frame())
sigma.e2.hat <- crossprod(fix.eff\$residual) / (fix.eff\$df.model)  #sigma_e^2
#Between effect model
cl3[[1]]<-as.name("bem")
between.eff <- eval(cl3,parent.frame())
RMSE2 <- crossprod(between.eff\$residual)/between.eff\$df2  #sigma_u^2
# T[i]
idx <- fix.eff\$id
Ti <- tapply(idx,idx, length)
g <- range(Ti)
Tbat <- length(Ti)/sum(1/Ti)
#default method for u2
sigma.u2.hat <- max(0,(RMSE2-sigma.e2.hat/Tbat))
##
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "id", "data"), names(mf), 0)
mf <- mf[c(1, m)]
mf[[1]]<-as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
X<-model.matrix(mt,mf);Y<-model.response(mf);id<-model.extract(mf,"id")
N <- dim(X)[1] #total obs
# use theta and partially dmean
Ti <- tapply(id,factor(id),length)[factor(id)]
theta <- 1- sqrt(sigma.e2.hat/(as.vector(Ti)*sigma.u2.hat+sigma.e2.hat))
#quasi-demean
xd<-dmean(X,id,theta=theta)
yd<-dmean(as.matrix(Y),id,theta=theta)
n<-dim(xd\$xm)[1]
Xbat <- xd\$xdm
Ybat <- yd\$xdm
fm <- lm.fit(Xbat,Ybat)
coe<-coef(fm) ;residual <- resid(fm)  ;fit <- fitted(fm)
sigma.e <- drop(crossprod(residual)/fm\$df.residual)
vcov<-sigma.e*solve(crossprod(Xbat))
std.Err<-matrix(sqrt(diag(vcov)),ncol=1)
rownames(std.Err) <- colnames(vcov)
#R2
R2.overall <- cor(X%*%coe,Y)^2
R2.between <-  cor(fix.eff\$xd\$xm %*% coe, fix.eff\$yd\$xm,use="pairwise.complete.obs")^2
R2.within <- cor(fix.eff\$xd\$xdm%*%coe,fix.eff\$yd\$xdm,use="pairwise.complete.obs")^2
##df
df2 <- fm\$df.residual
df1 <- N - df2 - 1
## wald test of the model The following code is wrong!!
#mdf <- ifelse(identical(g[1],g[2]),df1,df1-1)
#Wald <- F.model*mdf
g <- range(Ti)
#t-value and p
tval<-coe/std.Err
p.coef<-2*pnorm(abs(tval),lower=FALSE)
#result
result <- list(coef=coe,vcov=vcov,std.Err=std.Err,R2.overall=R2.overall,
R2.between=R2.between,R2.within=R2.within,g=g,call=cl,N=N,
n=n,tval=tval,p.coef=p.coef,sigma.e2.hat=sigma.e2.hat,
sigma.u2.hat=sigma.u2.hat,theta=theta)
class(result) <- "rem" ;return(result)
}
```

## Try the bpm package in your browser

Any scripts or data that you put into this service are public.

bpm documentation built on May 2, 2019, 5:46 p.m.