1
$\begingroup$

I have this simple example, where I try to get to the same factor loading as nlme::gls

Prepare data

#number of individuals
NN=2
#time
TT=100
set.seed(1234)
#predictors
x=matrix(rnorm(NN*TT,0,0.01),ncol=NN)
#correlated residuals
eps_ucor=matrix(rnorm(NN*TT,0,0.01),ncol=NN)
eps_cor=eps_ucor%*%chol(matrix(c(1,0.8,0.8,1),NN,NN))
#true factor loading
beta=1
#response
y=beta*x+eps_cor
#prepare data.frame; coerce by rows, i.e. y=c(y[1,1],y[1,2],y[2,1],y[2,2],...)
data=data.frame(t=sort(rep(1:TT,NN)),n=rep(1:NN,TT),y=c(t(y)),x=c(t(x)))

First, let's calculate the factor loading "manually" (I will follow https://en.wikipedia.org/wiki/Generalized_least_squares)

#covariance of residuals
cov_eps=cov(eps_cor)
#this is the matrix \Omega in the wikipedia article
#note that I assume that there is no correlation in time!
omega=diag(TT)%x%cov_eps
#calculate inverse
omega_inv=diag(TT)%x%solve(cov_eps)

We can now calculate the factor loading

beta_gls_1=(t(data$x)%*%omega_inv%*%data$x)^(-1)*t(data$x)%*%omega_inv%*%data$y
beta_gls_1
#         [,1]
#[1,] 1.017732

Finally, let's calculate the same with nlme::gls We have to prepare a vector with normalized variances

var_vec=diag(cov_eps)
var_vec=var_vec/var_vec[1]
names(var_vec)=1:2

Now let's perform the regression

res_gls=gls(model=y~x+0,
    data=data,
    correlation=corSymm(value=cov2cor(cov_eps)[2,1],form=~n|t,fixed=TRUE),
    weights=varIdent(value=numeric(0),form=~1|n,fixed=var_vec[-1]),method="ML")

The gls factor loading is

res_gls$coef
#       x 
#1.015082

As you can see, there is a difference between beta_gls_1 and res_gls$coef. I know, the difference is rather small - but from my experience it usually pays off to get to the bottom of such discrepancies. My guess is that I have misunderstood how to correctly initialize the fixed correlation and variance structure in gls. Any pointers to understand the difference highly appreciated!

BTW:

R version 3.0.2 (2013-09-25)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] nlme_3.1-113

loaded via a namespace (and not attached):
[1] grid_3.0.2      lattice_0.20-24 tools_3.0.2 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.