Hello, I would like to use mapply to avoid using a loop but for some reason, I can't seem to get it to work. I've included copies of my code below. The first set of code uses a loop (and it works fine), and the second set of code attempts to use mapply but I get a "subscript out of bounds" error. Any guidance would be greatly appreciated. Xj, Yj, and Wj are also lists, and s2, TAU, and GAMMA are scalars.
 
Thank You.
 
# THIS WORKS USING THE LOOP
for (j in 1:J) {
V.tilde.j <- solve((1/s2)*t(Xj[[j]])%*%Xj[[j]] + solve(TAU))
# Not singular case:
if(round(det(t(Xj[[j]])%*%Xj[[j]]),8)!=0) {
 Beta.hat.j <- solve(t(Xj[[j]])%*%Xj[[j]])%*%t(Xj[[j]])%*%Yj[[j]]
 V.j <- s2*solve(t(Xj[[j]])%*%Xj[[j]])
 Lambda.j <- solve(solve(V.j) + solve(TAU))%*%solve(V.j)
 Beta.tilde.j <- Lambda.j%*%Beta.hat.j + (diag(P) - Lambda.j)%*%Wj[[j]]%*%GAMMA
} 
# Singular case
else {
 Beta.tilde.j <- V.tilde.j%*%((1/s2)*t(Xj[[j]])%*%Yj[[j]] + solve(TAU)%*%Wj[[j]]%*%GAMMA)
}
BETA.Js[[j]] <- t(rmnorm(1, mean=as.vector(Beta.tilde.j), V.tilde.j))
}
 
 
 
# THIS DOESN'T WORK USING MAPPLY
update.betas <- function(s2,Xj,Yj,TAU,Wj,GAMMA) {
 V.tilde.j <- solve((1/s2)*t(Xj[[j]])%*%Xj[[j]] + solve(TAU))
 # Not singular case:
 if(round(det(t(Xj[[j]])%*%Xj[[j]]),8)!=0) {
  Beta.hat.j <- solve(t(Xj[[j]])%*%Xj[[j]])%*%t(Xj[[j]])%*%Yj[[j]]
  V.j <- s2*solve(t(Xj[[j]])%*%Xj[[j]])
  Lambda.j <- solve(solve(V.j) + solve(TAU))%*%solve(V.j)
  Beta.tilde.j <- Lambda.j%*%Beta.hat.j + (diag(P) - Lambda.j)%*%Wj[[j]]%*%GAMMA
 } 
 # Singular case
 else {
  Beta.tilde.j <- V.tilde.j%*%((1/s2)*t(Xj[[j]])%*%Yj[[j]] + solve(TAU)%*%Wj[[j]]%*%GAMMA)
 }
 BETA.Js[[j]] <- t(rmnorm(1, mean=as.vector(Beta.tilde.j), V.tilde.j))
 
 return(Beta.tilde.j)
}
BETA.Js <- mapply(update.betas,s2,Xj,Yj,TAU,Wj,GAMMA,SIMPLIFY=FALSE)
 


      
	[[alternative HTML version deleted]]


