fnc.form.Aj.theta5M <- (deriv.alphaj.s){ # fucntion create variance covariance matrix for lambda for a single school
vector_len <- (deriv.alphaj.s)
num_alpha <- (-1 + (1 + 4 * 2 * vector_len)) / 2
d <- (deriv.alphaj.s[1:num_alpha])
row_length <- (d)[1]
first_column <- 2
covs <- deriv.alphaj.s(row_length + 1):(deriv.alphaj.s)]
covs_count <- 1
for ( (2:row_length - 1)) {
for (column (first_column:row_length)) {
d[row,column] <- covs[covs_count]
d[column,] <- covs[covs_count]
covs_count <- covs_count + 1
}
first_column <- first_column + 1
}
(d)
}
### variables needeed in this function (lambda.tilde.mle, beta.mle, l.Omega.inv, )
VA_estimation <- (lambda.tilde.mle, beta.mle, l.Omega.inv, J, Y, Z, n, M){
l.lambdaj<-()
(j 1:J){
temp.lambdaj.tilde <- lambda.tilde.mle,j]
temp.lambdaj <- fnc.form.Aj.theta5M(temp.lambdaj.tilde)
l.lambdaj[j]] <- temp.lambdaj
}
#-----------------------------------------------------------------
# Omega.inv.mle
# beta.mle
# var.beta.mle<- solve(t(Z)%*%Omega.inv.mle%*%Z)
#
# # tstat.mle
# v.beta.mle<-diag(var.beta.mle); v.beta.mle
# #sd.beta.mle<-sqrt(v.beta.mle)
# #tstat.mle<-beta.mle/sd.beta.mle
#-----------------------------------------------------------------
# gamma.mj
lst.gamma<-()
count_begin <- 0
count_end <- 0
v.mle <- Y - Z %*% beta.mle
for (j 1:J){
n.j <- n[j]
count_begin <- count_end + 1
nrows <- n[j] * M
count_end <- count_end + nrows
omega.j.inv <- l.Omega.inv[j]]
a <- ((l.lambdaj[j]])) %x% (iota.nj(n.j))
gamma.j <- a %*% omega.j.inv %*% v.mle[count_begin:count_end]
lst.gamma[j]] <- gamma.j
count_begin <- count_end + 1
#browser()
}
<- ((, lst.gamma))
#browser()
()
#-----------------------------------------------------------------
# match university ID with gamma.mj's
# va.out<-cbind(idj,gamma)
#
# #-----------------------------------------------------------------
# va.out$ave_cent <- (va.out$gamma.raz+va.out$gamma.lec + va.out$gamma.eng + va.out$gamma.cit + va.out$gamma.wcom)/M
#
# #----- weights of the eta's from the standard deviations
#
# # turn into zeros the negative entries:
#
#
# l.lambdaj <- lapply(l.lambdaj, function(x) ifelse(x < 0, 0, x) )
# #dim(lambda.tilde)<-c((M*(M-1)),J)
#
# va.out$w_cent<-0
# for (j in 1:J)
# {
# #va.out$w_cent[j]<- as.matrix(gamma[j,])%*% solve((as.matrix(l.lambdaj[[j]]))^(0.5)) %*% (iota.nj(M)) # if singular, solve() does not work
# va.out$w_cent[j]<- as.matrix(gamma[j,])%*% ginv((as.matrix(l.lambdaj[[j]]))^(0.5)) %*% (iota.nj(M)) # if invertable, solve() and ginv() give the same exact result!!!!
# }
#
# #output <- va.out[order(va.out[,names(va.out)=="cent"]),]
#
# output <- va.out[with(va.out, order(va.out$w_cent, decreasing=TRUE)),]
# va.plot<- output
#
# # Correlation tables: among univariate VA's
#
# out.qr<-out.qr[,c(1,2,4)]; names(out.qr)[3]<-"va.qr"; names(out.qr)[2]<-"nj.qr"
# out.cr<-out.cr[,c(1,4)]; names(out.cr)[2]<-"va.cr"
# out.en<-out.en[,c(1,4)]; names(out.en)[2]<-"va.en"
# out.cc<-out.cc[,c(1,4)]; names(out.cc)[2]<-"va.cc"
# out.wrc<-out.wrc[,c(1,4)]; names(out.wrc)[2]<-"va.wrc"
#
#
# univar.va<-merge(out.qr, out.cr,by="univ_id")
# univar.va<-merge(univar.va, out.en,by="univ_id")
# univar.va<-merge(univar.va, out.cc,by="univ_id")
# univar.va<-merge(univar.va, out.wrc,by="univ_id")
#
# univar.va[1:5,]
# univar.cor<-cor(as.matrix(univar.va[,c(3:(M+2))]),method ="spearman")
# univar.cor<-round(univar.cor, 3)
#
# # Correlation tables: among multivariate and univariate VA's
#
# names(univar.va)[1]<-"idj"
# all.va<-merge(univar.va, va.out,by="idj")
# names(all.va)
# allva.cor<-cor(as.matrix(all.va[,c(3:7,10:16)]), method="spearman")
# allva.cor<-round(allva.cor,3)
#
# # correlation between modules
# module.cor<-cor(cbind(data.y$mod_razona_cuantitativo_punt,data.y$mod_lectura_critica, data.y$mod_ingles_punt,data.y$mod_comp_ciudadanas_punt, data.y$mod_comunica_escrita_punt))
# module.cor<-round(module.cor,3)
}
#-------------------------------------------------------------------------------
#gamma estimation for univariate-------------------------------------------------------------------------------
Univar_VA_estimation <- (x, y, n, x2 = ) {
# within estimator
# NOTE: For model 1 the intercept, and for models 2 & 3 both the intercept and the COMP variable turn to zeros in WZ (!)
Z.w <- x
WZ <- (Within2(Z.w,n))
WY <- (Within2(y,n)) # wy = wx*beta + wu
pi.w.hat<- ginv((WZ)) %*% (WZ,WY)
###pi.w.hat
# zeta.sq.hat
u.w.hat <- y- (Z.w %*% pi.w.hat); #is.matrix(u.w.hat)
Wu <- (Within2(u.w.hat, n))
N <- (Z.w)
J <- (n)
K <- (Z.w)
zeta.sq.hat <- (u.w.hat, Wu)/(N - J - K) ### sigma 2
###zeta.sq.hat; N; J; K
# between estimator
# NOTE: Should remove the COMP variable because one of the covariates and its COMP are identical in BZ (!)
# define Z.b
Z.b <- ((1,(n)), x) ### with intercept
BZ <- (Between2(Z.b, n))
BY <- (Between2(y, n))
pi.b.hat<- ginv((BZ)) %*% (BZ,BY)
# tau.sq.hat
u.b.hat<- BY - (BZ%*%pi.b.hat)
n_vector <- (n)
K <- (Z.b)
A <- (u.b.hat)
B <- zeta.sq.hat * (((BZ,((1/n_vector)%*%BZ))%*%((BZ)))) ###B
<- zeta.sq.hat * (1/n) ; ###C
tau.sq.hat <- (A+B-)/(J-K)
# GLS estimator of pi.hat (Z should include the intercept and the COMP (if acc. to the model))
# V.hat
tau <- (tau.sq.hat)
zeta <- (zeta.sq.hat)
# create a block diagonal matrix: V
listV<-()
(i (1:(n))){
temp <- J.nj(n[i]) * tau[1] + Diagonal(n[i]) * zeta[1]
listV[i]] <- (temp)
}
V.hat.inv<-(bdiag(listV))
# pi.gls.hat
if (!(x2)) {
Z <- ((1,(n)),x2)
}
{Z <- ((1,(n)),x)}
pi.var <- (((Z,V.hat.inv))%*%Z)
pi.sd <- ((pi.var))
pi.gls.hat <- pi.var %*% ((Z,V.hat.inv) %*% y)
# value added
u.hat<- (y-Z%*%pi.gls.hat)
N <- (u.hat)
K<- (Z)
tau.L.N <- tau[1]*(L.n(n_vector))
va <- tau.L.N%*%V.hat.inv%*%u.hat ###; va ; N;J;K
# link to university id and name:
(va)
}
va_multiple <- (x, y, n, x2 = ) {
df_va <- 0
for (i (1 : (y[1,]))) {
y_single <- (y,i])
if (!(x2)) {
va <- Univar_VA_estimation(x, y_single, n, x2)
}
{va <- Univar_VA_estimation(x = x, y = y_single, n = n)}
df_va <- (df_va, va)
}
(df_va)
}
