IGLS Algorithm

IGLS - 2l model

Step 1 - Estimate beta by OLS

zt<-matrix(1,nt)
vnaive<- diag(1,nt)   #iterate 1: sigma2=1
(beta<- solve(t(zt)%*%solve(vnaive)%*%zt)%*%t(zt)%*%solve(vnaive)%*%basef.f$score)

##          [,1]
## [1,] 248.4361

Step 2 - Calculate naive residual

res2<-basef.f$score-beta  

## Warning in basef.f$score - beta: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

Step 3 - Stack vectors from the matrix of cross product of residuals

y1<-tcrossprod(res2)
y<-as_tibble(as.vector(y1))

## Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
## This warning is displayed once per session.

Step 4 - Rewrite the covariance matrix as a regression equation

(t1<-as.vector(table(basef.f$id.school)))

## [1] 33 20

sum(t1)

## [1] 53

n1<-t1[1]
v11<-matrix(1,n1,n1)
v12<-diag(rep(1,n1))
v1<-v11+v12

n2<-t1[2]
v21<-matrix(1,n2,n2)
v22<-diag(rep(1,n2))
v2<-v21+v22

v<-adiag(v1,v2)
ver<-v%x%v # kronecker product
dim(ver)

## [1] 2809 2809

v<-as_tibble(as.vector(v))
dim(v)

## [1] 2809    1

z1<-mutate(v, 
           sigma_nu0 = if_else(value==0,0,1))
z1<-select(z1, sigma_nu0)
z1<-pull(z1)

z2<-mutate(v, 
           sigma_epsilon = if_else(value==2,1,0)) 
z2<-select(z2, sigma_epsilon)
z2<-pull(z2)

z3<-cbind(z1,z2)
dim(z3)

## [1] 2809    2

Step 5 - Estimate the residual variances

s_ver<-solve(ver)
theta1<- solve(t(z3)%*%s_ver%*%z3)#
theta2<-theta1%*%t(z3)%*%s_ver#
(theta<-theta2%*%as.vector(y1))

##         [,1]
## z1  227.8514
## z2 1823.5736

Step 7 - Estimate covariance matrix of residuals

(covf<- 2*theta1)

##              z1           z2
## z1  1.081686867 -0.001567251
## z2 -0.001567251  0.039215548