R/04_generalfit_FisherInformation.R

In SachaEpskamp/psychonetrics: Structural Equation Modeling and Confirmatory Network Analysis

numeric_FisherInformation <- function(model){
  model <- expectedmodel(model)
  # 2 * sum(model@sample@groups$nobs) * numDeriv::jacobian(psychonetrics_gradient,parVector(model),model=model)
  
  # Unit information instead:
  0.5 * numDeriv::jacobian(psychonetrics_gradient,parVector(model),model=model)
}

# General Fisher information former:
psychonetrics_FisherInformation <- function(model, analytic = TRUE){
  # Maybe do numeric?
  if (!analytic){
    return(numeric_FisherInformation(model))
  }

  # Prepare model:
  if (model@cpp){
    prep <- prepareModel_cpp(parVector(model), model) # <- upated!
  } else {
    prep <- prepareModel(parVector(model), model)  
  }
  
  
  # # estimator part Jacobian:
  # estimatorJacobian <- switch(
  #   model@estimator,
  #   "ML" = switch(model@distribution,
  #                 "Gaussian" = jacobian_gaussian_sigma
  #   )
  # )
  # estimatorPartJacobian <- estimatorJacobian(prep)
  
  # estimator part expected Hessian:
  if (model@cpp){
    
    estimatorHessian <- switch(
      model@estimator,
      "ML" = switch(model@distribution,
                    "Gaussian" = expected_hessian_Gaussian_cpp, # <- Updated!
                    "Ising" = expected_hessian_Ising
      ),
      "ULS" = switch(model@distribution,
                     "Gaussian" = expected_hessian_ULS_Gaussian_cpp # <- Updated!
      ),
      "WLS" = switch(model@distribution,
                     "Gaussian" = expected_hessian_ULS_Gaussian_cpp # <- Updated!
      ),
      "DWLS" = switch(model@distribution,
                      "Gaussian" = expected_hessian_ULS_Gaussian_cpp # <- Updated!
      ),
      "FIML" = switch(model@distribution,
                      "Gaussian" = expected_hessian_fiml_Gaussian_cppVersion # <- Updated!
      )
    )
    
  } else {
    estimatorHessian <- switch(
      model@estimator,
      "ML" = switch(model@distribution,
                    "Gaussian" = expected_hessian_Gaussian,
                    "Ising" = expected_hessian_Ising
      ),
      "ULS" = switch(model@distribution,
                     "Gaussian" = expected_hessian_ULS_Gaussian
      ),
      "WLS" = switch(model@distribution,
                     "Gaussian" = expected_hessian_ULS_Gaussian
      ),
      "DWLS" = switch(model@distribution,
                      "Gaussian" = expected_hessian_ULS_Gaussian
      ),
      "FIML" = switch(model@distribution,
                      "Gaussian" = expected_hessian_fiml_Gaussian
      )
    )
  }



  estimatorPart <- estimatorHessian(prep)
  
  # model part:
  if (model@cpp){
    modelJacobian <- switch(
      model@model,
      # "lnm" = d_phi_theta_lnm,
      # "ggm" = d_phi_theta_ggm,
      # "rnm" = d_phi_theta_rnm,
      # "gvar" = ifelse(model@rawts,d_phi_theta_gvar_rawts,d_phi_theta_gvar),
      "varcov" = d_phi_theta_varcov_cpp, # <- updated!
      "lvm" = d_phi_theta_lvm_cpp, # <- updated!
      "var1" = d_phi_theta_var1_cpp, # <- updated!
      # "panelvar1" = d_phi_theta_panelvar1,
      "dlvm1" = d_phi_theta_dlvm1_cpp, # <- updated!
      "tsdlvm1" = d_phi_theta_tsdlvm1_cpp, # <- updated!
      "meta_varcov" = d_phi_theta_meta_varcov_cpp, # <- updated!
      "Ising" = d_phi_theta_Ising_cpp, # <- updated!
      "ml_lvm" = d_phi_theta_ml_lvm_cpp # <- updated!
      # "cholesky" = d_phi_theta_cholesky
    )    
  } else {
    modelJacobian <- switch(
      model@model,
      # "lnm" = d_phi_theta_lnm,
      # "ggm" = d_phi_theta_ggm,
      # "rnm" = d_phi_theta_rnm,
      # "gvar" = ifelse(model@rawts,d_phi_theta_gvar_rawts,d_phi_theta_gvar),
      "varcov" = d_phi_theta_varcov,
      "lvm" = d_phi_theta_lvm,
      "var1" = d_phi_theta_var1,
      # "panelvar1" = d_phi_theta_panelvar1,
      "dlvm1" = d_phi_theta_dlvm1,
      "tsdlvm1" = d_phi_theta_tsdlvm1,
      "meta_varcov" = d_phi_theta_meta_varcov,
      "Ising" = d_phi_theta_Ising,
      "ml_lvm" = d_phi_theta_ml_lvm
      # "cholesky" = d_phi_theta_cholesky
    )
  }
  
  modelPart <- sparseordense(modelJacobian(prep))
  
  
  # Manual part:
  if (model@cpp){
    manualPart <- Mmatrix_cpp(model@parameters)
  } else {
    manualPart <- Mmatrix(model@parameters)  
  }
  # Compute fisher information and return:
  # Fisher <- 2 * prep$nTotal * t(manualPart) %*% t(modelPart) %*% estimatorPartHessian %*% modelPart %*% manualPart

  # Unit information instead:
  if (model@cpp){
   
    if (is(modelPart, "sparseMatrix")){
      Fisher <- FisherInformation_inner_cpp_DSS(as.matrix(estimatorPart), modelPart, manualPart)
    } else {
      Fisher <- FisherInformation_inner_cpp_DDS(as.matrix(estimatorPart), as.matrix(modelPart), manualPart)
    }
    
  } else {
    Fisher <- 0.5 * t(manualPart) %*% t(modelPart) %*% estimatorPart %*% modelPart %*% manualPart
  }
  
  
  as.matrix(Fisher)
}
# 
# 
# # Fisher information for LNM model
# Fisher_lnm <- function(model){
#   # Prepare
#   prep <- prepare_lnm(parVector(model), model)
# 
#   # d_phi_theta:
#   d_phi_theta <- d_phi_theta_lnm(prep)
#   
#   # Model matrix:
#   M <- Mmatrix(model@parameters)
#   
#   # Number of groups:
#   nGroups <- nrow(model@sample@groups)
#   
#   # Total number of parameters:
#   nTotal <- nrow(M)
#   nPar_perGroup <- nTotal / nGroups
#   
#   # number of constrained parameters:
#   nConstrained <- ncol(M)
#   
#   # Number of variables:
#   nVars <- nrow(model@sample@variables)
#   
#   # n means:
#   nMeans <- nVars
#   
#   # N var/covs:
#   nVarCov <- nVars*(nVars+1)/2
#   
#   # total per group:
#   nPerGroup <- nMeans + nVarCov
#   
#   # Create empty fisher information matrix:
#   I <- Matrix(0,nTotal,nTotal)
#   # For each group
#   for (g in 1:nGroups){
# 
#     # Indices of the mean part:
#     meanPart <- (g-1)*nPerGroup + seq_len(nMeans)
#     varPart <- (g-1)*nPerGroup + nMeans + seq_len(nVarCov)
#     groupPart <- (g-1)*nPar_perGroup + seq_len(nPar_perGroup)
#     kappa <- prep$groupModels[[g]]$kappa
#     
#     # Duplication mat:
#     D <- prep$groupModels[[g]]$D
# 
#     # Fill in
#     I[(g-1)*nPar_perGroup + seq_len(nPar_perGroup),(g-1)*nPar_perGroup + seq_len(nPar_perGroup)] <-
#       # (prep$nPerGroup[g] / prep$nTotal) *
#       ( prep$nTotal / prep$nPerGroup[g]) *
#       2*t(d_phi_theta[meanPart,groupPart]) %*% kappa %*% d_phi_theta[meanPart,groupPart] +
#       t(d_phi_theta[varPart,groupPart]) %*% t(D) %*% (kappa %x% kappa) %*% D %*% d_phi_theta[varPart,groupPart]
# 
# # 
# #     
# #     # Fill the relevant elements:
# #     # for (i in (g-1)*nPar_perGroup + seq_len(nPar_perGroup)){
# #     #   # Derivatives of sigma for i:
# #     #   dSigmai <- sparseMatrix(
# #     #     i = row(kappa)[lower.tri(kappa,diag=TRUE)],
# #     #     j = col(kappa)[lower.tri(kappa,diag=TRUE)],
# #     #     x = d_phi_theta[varPart, i], symmetric = TRUE)
# #     #   
# #     #   for (j in ((g-1)*nPar_perGroup + 1):i){
# #     #     # Derivatives of sigma for j:
# #     #     dSigmaj <- sparseMatrix(
# #     #       i = row(kappa)[lower.tri(kappa,diag=TRUE)],
# #     #       j = col(kappa)[lower.tri(kappa,diag=TRUE)],
# #     #       x = d_phi_theta[varPart, j], symmetric = TRUE)
# #     #     
# #     #     I[i,j] <- I[j,i] <- 2 * t(d_phi_theta[meanPart,i,drop=FALSE]) %*% kappa %*% d_phi_theta[meanPart,j,drop=FALSE] + 
# #     #        sum(diag(kappa %*% dSigmai %*% kappa %*% dSigmaj ))
# #     #   }
# #     # }
#   }
#   
# 
#   
#   # Return fisher information:
#   Fish <- t(M) %*% I %*% M
#   
#   # Return:
#   return(Fish)
# }