How to change the setting of delta matrix (for complementarity)

[yaoyao7031]

In the context of the eMDC model, I've noticed that the delta matrix for substitutability/complementarity remains consistent across all observations. However, due to the complexity of my data construction, I need the flexibility to modify this matrix between observations.

To illustrate, consider a simplified example. For the first observation, the matrix should be:
(0, 0, 0
d12, 0, 0
d13, d23, 0)
But for the second observation, I require a different matrix because the order of the alternatives are changed (although my actual data is more complicated than this)
(0, 0, 0
d13, 0, 0
d12, d23, 0)

I am aware that reorganizing the data could potentially align the matrices, but given the constraints of my research question, such reorganization is challenging. I'm wondering if the eMDC code can accommodate the adjustment of the delta matrix across different observations?

[dpalma]

Hi,

The delta matrix can be defined with different values for each observation in your sample.

Let's start with a situation with 3 alternatives, where everyone has the same complementarity/substitution (d) parameters. You would define the delta matrix as:

Code: Select all

delta <- matrix(c(  0,   0,   0, 
                  d21,   0,   0, 
		  d31, d32,   0), nrow=3, ncol=3, byrow=TRUE)

Where d21, d31, and d32 are the complementarity/substitution parameters for alternative pairs (2,1), (3,1), and (3,2).

Now let's imagine you have low and high income people in your sample, and you want them to have a different level of compl/subst for the (2,1) pair of alternatives. You would then have to define the delta "matrix" as a list of lists, as follows.

Code: Select all

delta <- list(
  list(                         0,   0, 0), 
  list(d21L*lowInc + d21H*highInc,   0, 0), 
  list(                       d31, d32, 0)
)

Where lowInc and highInc are dummy variables.

You can make the expressions inside these lists as complicated as you want. But I would recommend not including randomness inside delta (as in random parameters or additional error terms). Randomness inside delta might work, but we haven't tested it, so I would avoid it for now. But systematic variation of preferences (deterministic variability) is fine.

Best wishes
David

[yaoyao7031]

Hi David,

Thank you so much for your quick response. I really appreciate it.

I understand that systematic variation of preferences is allowed (say the complementarity/substitution patterns vary between different consumers based on their demographics). My question is a bit more complicated than that. Say there are 30 alternatives in total in the dataset, but each observation has a set of 3 alternatives to choose from. The baseline utility (beta b) and satiation parameters (gamma r) cannot be assigned uniquely to every single alternative.

To address this, I propose categorizing the 30 alternatives into three groups. Alternatives 1-10 are grouped as category 1, 11-20 as category 2, and 21-30 as category 3. This categorization simplifies the problem by allowing baseline utility (V) and satiation parameter (r) to be category-specific, rather than specific to each alternative. I still consider a budget allocation problem to 30 potential alternatives (although each observation only has 3 to allocate to), but baseline utility b and satiation parameter r are category-specific, not alternative-specific. I could run MDCEV model with some adjustments here.

Problem occurs when I try to run the eMDC model. The challenge arises in fitting it to the distinct characteristics of each observation. The conventional package assumes a universal delta matrix across all observations, with the possibility of demographic adjustments. However, my research context requires a unique delta matrix for each observation based on the category to which the alternatives belong.

For example, if Obs 1 faces alternatives 1, 11, 21 and Obs 2 faces alternatives 1, 11, 14, the delta matrix for Obs 1 is:
(0, 0, 0
d12, 0, 0
d13, d23, 0)
and that for Obs 2 is:
(0, 0, 0
d12, 0, 0
d12,0, 0)
This adjustment becomes essential due to the shared category between alternatives 11 and 14 in Observation 2.

That is why, I am wondering whether the package allows me to give a different delta matrix for each observation based on the set of alternatives/categories each observation faces. I hope my question makes sense. Thank you!

[dpalma]

Hi,

You can use the same trick for your situation, except that instead of using sociodemographics to determine the delta parameter to use in each case, you use the product category.

So let's take your situation. You have 30 possible products (products 1 to 30), but in each observation a decision maker only faces three alternatives: A, B, and C. Each alternative could be any of the 30 products. I would do the following:

Code: Select all

dBA <- d21*(A1*B2 + A2*B1) + d31*(A3*B1 + A1*B3) + d32*(A3*B2 + A2*B3)
dCA <- d21*(A1*C2 + A2*C1) + d31*(A3*C1 + A1*C3) + d32*(A3*C2 + A2*C3)
dCB <- d21*(B1*C2 + B2*C1) + d31*(B3*C1 + B1*C3) + d32*(B3*C2 + B2*C3)

delta <- matrix(c(  0,   0,   0, 
                  dBA,   0,   0, 
                  dCA, dCB,   0), nrow=3, ncol=3, byrow=TRUE)

Where A1 is a dummy variable taking value 1 if alternative A belongs to category 1. A2 is a dummy taking value 1 if alternative A belongs to category 2, and so on. You can create these variables in your database prior to estimation with something like:

Code: Select all

A1 <- prod_A %in%  1:10
A2 <- prod_A %in% 11:20
A3 <- prod_A %in% 21:30
B1 <- prod_B %in%  1:10
B2 <- prod_B %in% 11:20
B3 <- prod_B %in% 21:30
C1 <- prod_C %in%  1:10
C2 <- prod_C %in% 11:20
C3 <- prod_C %in% 21:30

Where prod_A is the product that alternative A is (so a number between 1 and 30).

I haven't double check the code, but the idea is there. I hope this helps.

Best wishes
David

[yaoyao7031]

Hi David,

Thank you very much for your reply. I tried your approach and it works!

I have another question related to this model. I tried to estimate a more complicated model (with EMDC model embedded in a latent class framework). I tried to use BFGS estimation. Everything worked fine in the estimation step and I received this information:

INFORMATION: Your model took more than 10 minutes to estimate, so it was saved to file XXX.rds before calculating its covariance matrix. If calculation of the covariance matrix fails or is stopped before finishing, you can load the model up to this point using apollo_loadModel. You may also want to inspect the approximate BHHH standard errors shown above to determine whether you wish to continue this process.

The BHHH s.e. looked fine here. However, after the numerical derivation of the covariance matrix, I saw that several s.e. (for the baseline utility parameters for one of the two segments) became NAN. I am wondering whether I should use the BHHH s.e. instead?

I also tried to use BHHH estimaton directly, but after about 10 iterations, I received the error message

Error in maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL, : NA in gradient

I am wondering what might be the problem here. I have 5 alternatives plus one outside good, so all the 5 baseline utilities should be identified. It also passed the "testing influence of parameters" stage. It means the model is not misspecified. Am I right?

[dpalma]

Hi,

Sorry for the delayed response. A few points before getting to your actual question.

• The fact that a model passes the "Testing the influence of parameters" stage does not imply that the model is correctly identified, as this is a necessary but not sufficient condition for identification.
• More generally, obtaining estimates is no proof of a model being identified. Many times you can get estimates without the model being identified.
• BHHH s.e. estimates are just an approximation of the real s.e.. Right now we do not know how accurate they are, so they are only there for reference, but you should not use them when reporting results in a publication, dissertation, consulting work, etc.
• The latest version of Apollo (v0.3.0) incorporates the new estimation algorithm BGW. I recommend using that over BFGS and BHHH.

Now, without much evidence, I believe yours is a problem of empirical lack of identification. In other words, you might be asking too much of your data. The eMDC model is already quite demanding in terms of identification: you have base utilities, satiation, and s.d. parameters to estimate. I have never estimated this model inside an ICLV myself. So while the code does support mixing, I am not sure how it behaves in terms of identification. My guess is it would at least require a good amount of observations.

Furthermore, where are you including the latent variable? I would recommend including mixing only in the base utility, not in the satiation parameters.

Best wishes
David

[yaoyao7031]

Hi David,

Thank you very much for your reply! I have been doing some debugging. I now have a few follow-up questions/clarifications:

1. In my panel data, I have about 500 individuals and about 25000 observations in total. Instead of using random coefficients, I use a latent class framework similar to Kamakura and Russell (1989) to capture heterogeneity, so I do not have s.d. parameters to estimate. Instead, if I try a two-segment model, I have two sets of baseline utility and satiation intercepts (alternative specific) to estimate. I also added the same set of covariates to both baseline and satiation. I think given the sample size, it should not have much identification issue to estimate about 60 parameters?

2. Following Kamakura and Russell's framework, I use some observed demographic variables to model the probability of each individual belonging to which segment. Thus, each segment has a different set of parameters for satiation and baseline intercepts and covariates. Is this what you asked about where I included latent variables?

3. When I tried different starting values using "apollo_searchStart", most of the starting values were rejected. I found the problem comes from the apollo_emdc1 function:

Code: Select all

G <- mapply(function(pi, Ei) phi0x0 * pi - Ei, cost,  E, SIMPLIFY = FALSE)
W <- mapply(function(vi, xi, gi, Gi) (vi - log(xi/gi + 1) - log(Gi)), V, X, gamma, G, SIMPLIFY = FALSE)

I believe W is the same as Equation (10) in Palma and Hess (2022). However, G, the component within the log, could be negative when the delta_kl parameters take certain values. In my case, phi0x0 * pi is relatively small. If delta_kl is a bit larger, G becomes negative and cannot take log. That is why when I estimate the model, it shows a lot of warnings saying

In log(Gi) : NaNs produced
In log(Jdet) : NaNs produced

I am wondering whether this was the reason that the estimation does not converge, and some times to standard errors are reported. Is there any way to make sure that the G would be positive when estimating using the package?

Thank you and look forward to your insights

[dpalma]

Hi,

Sorry for the delayed response.

Kamakura & Russell (1989)’s approach is what most people in choice modelling would call a latent class approach. So using their notation you have the f_ni probability, which is the class allocation probability, so the probability of individual n belonging to class i.

In Kamakura & Russell, they only use constants (lambda) to define the probability of each individual belonging to a different class, so all individuals have the same probability of belonging to each class. But you could further parametrise lambda, so that it depends on the characteristic of the individual. This way, f_ni could be different for each individual n. People usually put sociodemographics as explanatory variables of f_ni, so you can later interpret each segment as that of the younger or older people, or the richer people, etc. Of course, the allocation is still random, so this means that richer or older people have more probability of belonging to a segment, but the allocation is not deterministic (as in clustering). And just as you can include sociodemographics to explain class allocation (f_ni) you can also use latent variables.

In terms of your specific question:

1. You seem to have enough observations to estimate a model with two (or maybe more) segments or classes. Keep in mind, though, that segmented or latent class models are not easy to estimate. The solution you get can be highly dependent on your starting point. So, I would recommend testing several starting points to make sure you get to a good solution.

2. As I mentioned before, you can make f_ni a function of characteristics of the individuals, including their sociodemographics, or latent variables that you may have calculated. Implementation in apollo gets a bit more complicated if you use latent variables in f_ni, but it can be done.

3. This is a known difficulty when estimating the model. For theoretical reasons, we want to allow marginal utilities to become negative under certain circunstances (due to the substitution effects), so there is no way to force Gi to be positive. That said, the optimisation algorithm should naturally avoid those areas, as when it runs into a solution that is not defined (where Gi < 0), it will test other values. So the warnings are not important in this case. But again, you might want to try different starting values. But throughout all the starting values you try, I would recommend keeping deltas always starting as zeros, to minimise this kind of problems.

Hope this helps clarify things a bit.

Best wishes
David