11# # File Name: invariance.alignment.R
2- # # File Version: 3.972
2+ # # File Version: 4.017
33
44
55invariance.alignment <- function ( lambda , nu , wgt = NULL ,
66 align.scale = c(1 ,1 ), align.pow = c(.5 ,.5 ), eps = 1e-3 ,
77 psi0.init = NULL , alpha0.init = NULL , center = FALSE , optimizer = " optim" ,
88 fixed = NULL , meth = 1 , vcov = NULL , eps_grid = seq(0 ,- 10 , by = - .5 ),
9- num_deriv = FALSE , ... )
9+ num_deriv = FALSE , le = FALSE , ... )
1010{
1111 CALL <- match.call()
1212 s1 <- Sys.time()
@@ -90,7 +90,8 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
9090 }
9191
9292 # -- define optimization functions
93- ia_fct_optim <- function (x , lambda , nu , overparam , eps , meth_ ){
93+ ia_fct_optim <- function (x , lambda , nu , overparam , eps , meth_ ,
94+ item_wise = FALSE ){
9495 res <- invariance_alignment_define_parameters(x = x , ind_alpha = ind_alpha ,
9596 ind_psi = ind_psi , reparam = reparam , meth = meth_ )
9697 alpha0 <- res $ alpha0
@@ -100,7 +101,11 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
100101 wgt = wgt , align_scale = align.scale ,
101102 align_pow = align.pow , eps = eps , wgt_combi = wgt_combi , type = type ,
102103 reparam = FALSE , meth = meth_ )
103- val <- val $ fopt
104+ if (item_wise ){
105+ val <- val $ fopt_item
106+ } else {
107+ val <- val $ fopt
108+ }
104109 if (overparam | meth == 0 ){
105110 G <- nrow(lambda )
106111 fac <- sum(wgt [,1 ]) / 1000
@@ -201,25 +206,41 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
201206 }
202207
203208 # -- standard error computation (if requested)
209+ TEjkbc <- LEjkbc <- TEjk <- LEjk <- TE <- TEbc <- LEbc <- LE <- SE <- NULL
210+ vcov0 <- vcov
211+ V_TEjkbc <- V_TEjk <- V_LEjk <- V_TEbc <- V_TE <- V_LEbc <- V_LE <- NULL
212+ V_LEjkbc <- NULL
213+
204214 if (! is.null(vcov )){
205215
206216 names(par ) <- c( paste0(' alpha' ,2L : G ), paste0(' psi' ,2L : G ) )
207217 NP <- length(par )
208218
209219 # - gradient computation
210- ia_grad_optim_num <- function (x , lambda , nu , overparam , eps , meth = 1 , h = 1e-4 ){
220+ ia_grad_optim_num <- function (x , lambda , nu , overparam , eps , meth = 1 ,
221+ h = 1e-4 , item_wise = FALSE ){
211222 NP <- length(x )
212223 par <- x
213- grad <- rep(0 ,NP )
214- names(grad ) <- names(x )
224+ if (! item_wise ){
225+ grad <- rep(0 ,NP )
226+ names(grad ) <- names(x )
227+ } else {
228+ grad <- matrix (0 ,nrow = NP , ncol = I )
229+ rownames(grad ) <- names(x )
230+ }
215231 args <- list (x = par , lambda = lambda , nu = nu , overparam = overparam , eps = eps ,
216- meth_ = meth )
232+ meth_ = meth , item_wise = item_wise )
217233 for (pp in 1L : NP ){
218234 args $ x <- mgsem_add_increment(x = par , h = h , i1 = pp )
219235 f1 <- do.call(what = ia_fct_optim , args = args )
220236 args $ x <- mgsem_add_increment(x = par , h = - h , i1 = pp )
221237 f2 <- do.call(what = ia_fct_optim , args = args )
222- grad [pp ] <- (f1 - f2 )/ (2 * h )
238+ der <- (f1 - f2 )/ (2 * h )
239+ if (! item_wise ){
240+ grad [pp ] <- der
241+ } else {
242+ grad [pp ,] <- der
243+ }
223244 }
224245 return (grad )
225246 }
@@ -265,11 +286,97 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
265286 }
266287 }
267288
268- A <- MASS :: ginv(hess_par ) %*% hess_theta
269- vcov <- A %*% vcov %*% t(A )
289+ H1 <- MASS :: ginv(hess_par )
290+ A <- H1 %*% hess_theta
291+ vcov <- A %*% vcov0 %*% t(A )
270292 rownames(vcov ) <- colnames(vcov ) <- names(par )
293+ SE <- sqrt_diag(x = vcov )
294+
295+ # --- linking error based on M-estimation
296+ if (le ){
297+ args1 <- list (x = par , lambda = lambda , nu = nu , overparam = overparam , eps = eps ,
298+ meth = meth , item_wise = TRUE )
299+
300+ # item-wise gradient
301+ grad_item <- do.call( what = ia_grad_optim_num , args = args1 )
302+ # variance matrix (meat matrix) for linking error estimation
303+ M <- 0 * vcov
304+ Mbc <- M
305+ for (ii in 1L : I ){
306+ gii <- grad_item [,ii ]
307+ Mii <- outer(gii ,gii )
308+ ind <- seq(ii ,2 * I * G ,I )
309+ vcov0_ii <- vcov0 [ind , ind ]
310+ hess_theta_ii <- hess_theta [,ind ]
311+ Mbcii <- hess_theta_ii %*% vcov0_ii %*% t(hess_theta_ii )
312+ M <- M + Mii
313+ Mbc <- Mbc + Mbcii
314+ }
315+ V_LE <- I / (I - 1 ) * H1 %*% M %*% t(H1 )
316+ V_LEbc <- I / (I - 1 ) * H1 %*% (M - Mbc ) %*% t(H1 )
317+ LE <- sqrt_diag(x = V_LE , names = names(par ))
318+ LEbc <- sqrt_diag(x = V_LEbc , names = names(par ))
271319
272- }
320+ V_TE <- vcov + V_LE
321+ V_TEbc <- vcov + V_LEbc
322+ TE <- sqrt(SE ^ 2 + LE ^ 2 )
323+ TEbc <- sqrt(SE ^ 2 + LEbc ^ 2 )
324+
325+ # --- jackknife linking error
326+
327+ hess_par_item <- list ()
328+ for (ii in 1L : I ){
329+ hess_par_item [[ii ]] <- 0 * hess_par
330+ }
331+
332+ H0 <- 0 * hess_par
333+ pp <- 1
334+ for (pp in 1L : NP ){
335+ args1 $ x <- mgsem_add_increment(x = par , h = h , i1 = pp )
336+ f1 <- do.call( what = ia_grad_optim_num , args = args1 )
337+ args1 $ x <- mgsem_add_increment(x = par , h = - h , i1 = pp )
338+ f2 <- do.call( what = ia_grad_optim_num , args = args1 )
339+ der <- (f1 - f2 )/ (2 * h )
340+ for (ii in 1L : I ){
341+ hess_par_item [[ii ]][,pp ] <- der [,ii ]
342+ }
343+ }
344+ for (ii in 1L : I ){
345+ H0 <- H0 + hess_par_item [[ii ]]
346+ }
347+ # hess_par=H0
348+
349+ # jk le
350+ jkfac <- (I - 1 )/ I
351+ estdiff_jk <- matrix (NA , nrow = I , ncol = NP )
352+ V_LEjkbc <- 0 * V_LE
353+
354+ for (ii in 1L : I ){
355+ cii <- grad_item [,ii ]
356+ Aii <- hess_par - hess_par_item [[ii ]]
357+ A1ii <- MASS :: ginv(Aii )
358+ parii <- A1ii %*% cii
359+ estdiff_jk [ii ,] <- parii [,1 ]
360+ # correction: see LESL paper
361+ U <- A
362+ ind <- seq(ii , 2 * I * G , I )
363+ HTii <- hess_theta
364+ HTii [,ind ] <- 0
365+ Ui <- A1ii %*% HTii
366+ Wi <- Ui - U
367+ V_LEjkbc <- V_LEjkbc + jkfac * ( Wi %*% vcov0 %*% t(Wi ) )
368+ }
369+ V_LEjk <- jkfac * crossprod(estdiff_jk )
370+ LEjk <- sqrt_diag(x = V_LEjk , names = names(par ))
371+ V_LEjkbc <- V_LEjk - V_LEjkbc
372+ LEjkbc <- sqrt_diag(x = V_LEjkbc , names = names(par ))
373+
374+ TEjk <- sqrt( SE ^ 2 + LEjk ^ 2 )
375+ TEjkbc <- sqrt( SE ^ 2 + LEjkbc ^ 2 )
376+
377+ } # end le
378+
379+ } # end vcov
273380
274381 # center parameters
275382 res <- invariance_alignment_center_parameters(alpha0 = alpha0 , psi0 = psi0 ,
@@ -325,7 +432,10 @@ invariance.alignment <- function( lambda, nu, wgt=NULL,
325432 lambda.resid = lambda.resid , nu.aligned = nu.aligned , lambda = lambda0 ,
326433 nu = nu0 , nu.resid = nu.resid , fopt = fopt , align.scale = align.scale ,
327434 align.pow = align.pow0 , res_optim = res_optim , eps = eps , wgt = wgt ,
328- miss_items = missM , numb_items = numb_items , vcov = vcov ,
435+ miss_items = missM , numb_items = numb_items , vcov = vcov , V_LE = V_LE ,
436+ V_LEbc = V_LEbc , V_TE = V_TE , V_LEjk = V_LEjk , V_LEjkbc = V_LEjkbc ,
437+ SE = SE , LE = LE , LEbc = LEbc , LEjk = LEjk , LEjkbc = LEjkbc ,
438+ TE = TE , TEbc = TEbc , TEjk = TEjk , TEjkbc = TEjkbc ,
329439 fct_optim_inits = fct_optim_inits , fixed = fixed , meth = meth0 ,
330440 meth_internal = meth , s1 = s1 , s2 = s2 , time_diff = time_diff , CALL = CALL )
331441 class(res ) <- ' invariance.alignment'
0 commit comments