library(OpenMx) mxOption(key='Number of Threads', value=parallel::detectCores() - 1) # arrow = 1 from -> to # arrow = 2 from <-> to # free = FALSE -> fixed path (parameter values are ex ante fix) # free = TRUE is a free path (parameter values are estimated) # from: sources of the new paths # values: Starting values of parameters # connect: type of source to sink connection # labels: name of the paths # values is the starting value of a path # reasonable values for faster convergence # first try to replicate model in # "Linear dynamic panel-data estimation using maximum likelihood and structural equation modeling" # yit = λyit−1 + x∗itβ + w∗iδ + αi + ξt + υit ∀i, t # E(υit|yt−1,i , xt,i, wi, αi) = 0 ∀i, library(haven) wages <- haven::read_dta("https://www3.nd.edu/~rwilliam/statafiles/wageswide.dta") #View(wages) head(wages) summary(wages) ## 595 households who reported a non-zero wage (measured for 7 years: 1976 - 1982) # y = wks: number of weeks employed in each year # x = union (dummy): 1 if wage set by union contract # w = lwage: ln of wage in each year # z = ED: years of education in 1976 data <- mxData(wages, type = "raw") #mxTypes() latents <- "alpha" manifests <- c("wks1", "union1", "lwage1", "wks2", "union2", "lwage2", "wks3", "union3", "lwage3", "wks4", "union4", "lwage4", "wks5", "union5", "lwage5", "wks6", "union6", "lwage6", "wks7", "union7", "lwage7", "ed") cov_data <- cov(wages[,manifests]) #mxTypes() model1 <- mxModel("Test1", type = "RAM", # Manifest variables manifestVars = c("wks1", "union1", "lwage1", "wks2", "union2", "lwage2", "wks3", "union3", "lwage3", "wks4", "union4", "lwage4", "wks5", "union5", "lwage5", "wks6", "union6", "lwage6", "wks7", "union7", "lwage7", "ed"), # Latent variable latentVars = c("alpha"), #mxPath(from = "alpha", to = manifests, # free = FALSE, values = c(rep(1, 21), 0)), # Correlation between alpha and all exogenous variables, except "ed" mxPath(from = c("alpha"), to = c("wks1", "union1", "lwage1", "union2", "lwage2", "union3", "lwage3", "union4", "lwage4", "union5", "lwage5", "union6", "lwage6", "union7", "lwage7", "ed"), arrows = 2, connect = "unique.pairs", free = FALSE, values = c(rep(1, 15), 0)), # Correlation between alpha and endogenous variables, fixed effects mxPath(from = c("alpha"), to = c("wks2", "wks3", "wks4", "wks5", "wks6", "wks7"), arrows = 1, connect = "unique.pairs", free = FALSE, values = c(rep(1, 6))), # Correlation between exogenous and endogenous variables mxPath(from = c("ed", "wks1", "lwage1", "union1"), to = c("wks2"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("ed", "wks2", "lwage2", "union2"), to = c("wks3"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("ed", "wks3", "lwage3", "union3"), to = c("wks4"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("ed", "wks4", "lwage4", "union4"), to = c("wks5"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("ed", "wks5", "lwage5", "union5"), to = c("wks6"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("ed", "wks6", "lwage6", "union6"), to = c("wks7"), arrows = 1, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4)), labels = c("ed", "lambda", "beta2", "beta3")), mxPath(from = c("wks2"), to = c("union3", "union4", "union5", "union6"), arrows = 2, connect = "unique.pairs", free = TRUE, values = c(rep(1, 4))), mxPath(from = c("wks3"), to = c("union4", "union5", "union6"), arrows = 2, connect = "unique.pairs", free = TRUE, values = c(rep(1, 3))), mxPath(from = c("wks4"), to = c("union5", "union6"), arrows = 2, connect = "unique.pairs", free = TRUE, values = c(rep(1, 2))), mxPath(from = c("wks5"), to = c("union6"), arrows = 2, connect = "unique.pairs", free = TRUE, values = 1), # Residual variances mxPath(from = c("wks2", "wks3", "wks4", "wks5", "wks6", "wks7"), arrows = 2, connect = "single", free = TRUE, values = c(rep(1, 5)), labels = rep("residual", 6)), mxPath(from = "one", to = manifests, arrows = 1, free = TRUE, values = 1), #Data mxData(wages, type = "raw", numObs = 595)) model1 mxRun(model1)