The R package factormodel provides functions to estimate a factor model using either discrete or continuous proxy variables. Such model is useful when proxy variables include measurement errors.
When proxy variables are discrete, you can use ‘dproxyme’ function. The function estimates a finite-mixture model using an EM algorithm (Dempster, Laird, Rubin, 1977).
A function ‘dproxyme’ returns a list of estimated measurement (stochastic) matrices and a type probability matrix from discrete proxy variables. The ij-th element in a measurement matrix is the conditional probability of observing j-th (largest) proxy response value conditional on that the latent type is i. The type probability matrix is of size N (num of obs) by sbar (num of type). The ij-th element of the type probability is the probability of observation i to belong to the type j. For further explanation on identification of measurement stochastic matrices, see Hu(2008) and Hu(2017).
When proxy variables are continuous, you can use ‘cproxyme’ function. The function estimates a linear factor model assuming a continuous latent variable.
A function ‘cproxyme’ returns a list of linear factor model coefficients and the variance of measurement errors in each proxy variable. For further explanation on identification of linear factor model, see Cunha, Heckman, Schennach (2010).
You can install a package factormodel using either CRAN or github.
or
library(factormodel)
library(nnet)
library(pracma)
library(stats)
library(utils)
# DGP
# set parameters
nsam <- 5000
M1 <- rbind(c(0.8,0.1,0.1),c(0.1,0.2,0.7))
M2 <- rbind(c(0.7,0.2,0.1),c(0.2,0.2,0.6))
M3 <- rbind(c(0.9,0.05,0.05),c(0.1,0.1,0.8))
CM1 <- t(apply(M1,1,cumsum))
CM2 <- t(apply(M2,1,cumsum))
CM3 <- t(apply(M3,1,cumsum))
# 40% of sample is type 1, 60% is type 2
truetype <- as.integer(runif(nsam)<=0.4) +1
# generate fake data
dat <- data.frame(msr1=rep(NA,nsam),msr2=rep(NA,nsam),msr3=rep(NA,nsam))
for (k in 1:nsam){
dat$msr1[k] <- which(runif(1)<=CM1[truetype[k],])[1]
dat$msr2[k] <- which(runif(1)<=CM2[truetype[k],])[1]
dat$msr3[k] <- which(runif(1)<=CM3[truetype[k],])[1]
}
# estimate using dproxyme
oout <- dproxyme(dat=dat,sbar=2,initvar=1,initvec=NULL,seed=210313,tol=0.005,maxiter=200,miniter=10,minobs=100,maxiter2=1000,trace=FALSE)
# check whether the estimated measurement stochastic matrices are same with the true # measurement stochastic matrices
print(oout$M_param)
#> [[1]]
#> 2 3
#> [1,] 0.78073858 0.1129363 0.1063251
#> [2,] 0.09111655 0.2031201 0.7057633
#>
#> [[2]]
#> 2 3
#> [1,] 0.6987843 0.2011567 0.1000589
#> [2,] 0.2053359 0.1939312 0.6007329
#>
#> [[3]]
#> 2 3
#> [1,] 0.89619844 0.0485488 0.05525276
#> [2,] 0.08761801 0.1015417 0.81084033
# check type probability
print(head(oout$typeprob))
#> [,1] [,2]
#> [1,] 0.9932423 0.006757696
#> [2,] 0.9979306 0.002069438
#> [3,] 0.1077760 0.892223987
#> [4,] 0.9050997 0.094900295
#> [5,] 0.9932423 0.006757696
#> [6,] 0.4947413 0.505258684Below example shows how to use ‘cproxyme’ function to estimate a linear factor model. The code first simulates fake data using a data generating process provided below and then estimates the parameters using ‘cproxyme’ function.
library(factormodel)
library(stats)
library(utils)
library(gtools)
#>
#> Attaching package: 'gtools'
#> The following object is masked from 'package:pracma':
#>
#> logit
set.seed(seed=210315)
# DGP
# set parameters
nsam <- 5000 # number of observations
np <- 3 # number of proxies
true_mtheta <- 2
true_vartheta <- 1.5
true_theta <- rnorm(nsam, mean=true_mtheta, sd=sqrt(true_vartheta))
# first proxy variable is an anchoring variable
true_alpha0 <- c(0,2,5)
true_alpha1 <- c(1,0.5,2)
true_varnu <- c(0.5,2,1)
# simulate fake data
dat <- matrix(NA,nrow=nsam,ncol=np)
for (k in 1:np){
dat[,k] <- true_alpha0[k] + true_alpha1[k]*true_theta + rnorm(nsam,mean=0,sd=sqrt(true_varnu[k]))
}
# estimate parameters using cproxyme
oout <- cproxyme(dat=dat,anchor=1)
# print estimated parameters
print(oout$alpha0)
#> [1] 0.000000 2.032455 5.086428
print(oout$alpha1)
#> [1] 1.0000000 0.4913708 1.9630702
print(oout$varnu)
#> [1] 0.4827664 2.0430161 1.0605388
print(oout$mtheta)
#> [1] 1.990616
print(oout$vartheta)
#> [1] 1.586096This vignette showed how to use functions in `factormodel’ R package.