Marginsplot

[Ansuyah Magan]

Hi
I am having difficulty with a margins plot in Stata 14.1. I measured bone size for 3 bones at different ages. My commands were:

xtmixed size age bone bone#c.age
margins bone, dydx (age)
quietly margins bone, at(age=(0(1)10))

All the above worked.
I got an error message with the following command:

marginsplot, xdimension(at(yplage)) noci recast(line)

Error message: invalid at() dimension information; using variable bone as a factor variable and a regular variable is not supported


I understand that this is also a linear prediction. What would be the steps to fit a quadratic to the data?

With thanks
Ansuyah

[Raffaele Grotti]

Hi Ansuyah,
try to get rid of the at() within the xdimension() option.

Raffaele

[Clyde Schechter]

I think you just misspecified the argument in -xdimension()-. I assume you want the horizontal axis dimension to be age, so just write

marginsplot, xdimension(age) noci recast(line)

Actually, I don't think you even need to specify -xdimension()- here. I'm pretty sure that the default x-axis variable here will be age anyway.

As for fitting a model that is quadratic in age:

mixed size i.bone##c.age##c.age

Note: Unless you are using an old version of Stata, the command's correct name is now -mixed-, without the -xt- prefix. -xtmixed- still works for now, but may become unsupported in future versions of Stata. So best to get in the habit of calling it -mixed- now. Also, you are using the -mixed- command but you haven't specified any random effects. So this is pointless, you may as well just use -regress-.

[Marcos Almeida]

The term yplage is apparently a mistype.

[Ansuyah Magan]

Thank you all for your suggestions.
At some ages I have more than one observation. I managed to produce graph when I simplified this to a single observation at a time point. How do I go about retaining these additional observations? I understand that I am supposed to use the "weight" option but do not know how to do this.
"yplage" was a typo.
"xdimension" needed to be specified whilst removing "at()" still produced graph
I used xtmixed rather than regress because I adjusted for the random effects of subject - I didn't include it here for the sake of simplicity
"mixed" vs "xtmixed" yields the same coefficients for the 3 bones but the residuals for the random effects differ.

Regards
Ansuyah

[Marcos Almeida]

"mixed" vs "xtmixed" yields the same coefficients for the 3 bones but the residuals for the random effects differ.

You are probably mistaking standard deviation (the default with xtmixed) with variance (the default with mixed). If you dutifully demand the variance option with xtmixed, you may get similar results, as shown below:

. webuse pig
(Longitudinal analysis of pig weights)

.  mixed weight week || id:

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -1014.9268  
Iteration 1:   log likelihood = -1014.9268  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        432
Group variable: id                              Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

                                                Wald chi2(1)      =   25337.49
Log likelihood = -1014.9268                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |   6.209896   .0390124   159.18   0.000     6.133433    6.286359
       _cons |   19.35561   .5974059    32.40   0.000     18.18472    20.52651
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   14.81751   3.124226      9.801716    22.40002
-----------------------------+------------------------------------------------
               var(Residual) |   4.383264   .3163348      3.805112     5.04926
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 472.65        Prob >= chibar2 = 0.0000

. predict mymixresid*, reffects

.  xtmixed weight week || id:

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -1014.9268  
Iteration 1:   log likelihood = -1014.9268  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        432
Group variable: id                              Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

                                                Wald chi2(1)      =   25337.49
Log likelihood = -1014.9268                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |   6.209896   .0390124   159.18   0.000     6.133433    6.286359
       _cons |   19.35561   .5974059    32.40   0.000     18.18472    20.52651
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                   sd(_cons) |   3.849352   .4058119      3.130769    4.732866
-----------------------------+------------------------------------------------
                sd(Residual) |   2.093625   .0755472       1.95067    2.247056
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 472.65        Prob >= chibar2 = 0.0000

. predict myxtmixresid*, reffects

.  xtmixed weight week || id:, var

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -1014.9268  
Iteration 1:   log likelihood = -1014.9268  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        432
Group variable: id                              Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

                                                Wald chi2(1)      =   25337.49
Log likelihood = -1014.9268                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |   6.209896   .0390124   159.18   0.000     6.133433    6.286359
       _cons |   19.35561   .5974059    32.40   0.000     18.18472    20.52651
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   14.81751   3.124226      9.801716    22.40002
-----------------------------+------------------------------------------------
               var(Residual) |   4.383264   .3163348      3.805112     5.04926
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 472.65        Prob >= chibar2 = 0.0000

. list my* in 1/10

     +-----------------------+
     | mymixre~1   myxtmix~1 |
     |-----------------------|
  1. | -1.683105   -1.683105 |
  2. | -1.683105   -1.683105 |
  3. | -1.683105   -1.683105 |
  4. | -1.683105   -1.683105 |
  5. | -1.683105   -1.683105 |
     |-----------------------|
  6. | -1.683105   -1.683105 |
  7. | -1.683105   -1.683105 |
  8. | -1.683105   -1.683105 |
  9. | -1.683105   -1.683105 |
 10. |  .8987018    .8987018 |
     +-----------------------+

[Ansuyah Magan]

Hi Marcos
Your explanation was most helpful. Thank you.
Regards
Ansuyah

[Ansuyah Magan]

Please could you advise on how to correct or adjust for more than 1 observation at a particular time point.

Thanks
Ansuyah

[Katherine Adams]

Hello,

I have a question about the marginsplot command.

Here is the example of my data; the data in the example is sorted by location, so, in fact, it is related to only one household with location id 600001 (the original data is for many households over 2017-2018).

* Example generated by -dataex-. To install: ssc install dataex
clear
input long location float(lconsum tp) byte randomgrp float(calday treatgr treat_numb_of_days)
600001  4.342219 0 0 20820 0 -403
600001  4.396476 0 0 20821 0 -402
600001 4.4473995 0 0 20822 0 -401
600001 4.4349075 0 0 20823 0 -400
600001 4.3400753 0 0 20824 0 -399
600001  3.974716 0 0 20825 0 -398
600001 4.2170517 0 0 20826 0 -397
600001 4.4074755 0 0 20827 0 -396
600001 4.2367565 0 0 20828 0 -395
600001 4.3245976 0 0 20829 0 -394
600001 4.2221044 0 0 20830 0 -393
600001 4.4336513 0 0 20831 0 -392
600001  4.122668 0 0 20832 0 -391
600001  4.443582 0 0 20833 0 -390
600001 4.1282955 0 0 20834 0 -389
600001  4.153299 0 0 20835 0 -388
600001  4.019543 0 0 20836 0 -387
600001 3.8549176 0 0 20837 0 -386
600001  3.745078 0 0 20838 0 -385
600001  3.776974 0 0 20839 0 -384
end
format %td calday

location; household’s location id
lconsum; log of energy consumption
tp; post-treatment variable; gen tp = (calday >= td(08feb2018))
randomgr; one of three treatment groups (can be 1,2,3, as well as 0 if it is a control group)
calday; day and year 01jan2017
treatgr; treatment indicator;
gen treatgr = randomgrp if calday >= td(08feb2018)
replace treat = 0 if treat ==.
treat_numb_of_days; the number of days before/after the treatment date
gen treat_numb_of_days = calday-td(08feb2018)


I ran an event-study regression (as I understand it):
areg lconsum treatgr tp i.treatgr#c.treat_numb_of_days, absorb(location) vce(cluster location)

Now, I need to plot an event-study figure with the interaction (‘i.treat#c.treat_numb_of_day’) coefficient on the y-axis and the number of days before/after the treatment on the x-axis (‘treat_numb_of_days’):
margins, at(treat_numb_of_days=(-403(86)285)) over(treatgr)
marginsplot
[the value of ‘treat_numb_of_days’ is -403, -402, …, 284, 285 in the original data; I chose the increment of 86 to simplify]


However, running the command ‘marginsplot’ gives the error similar to that Ansuyah Magan had:
invalid at() dimension information; using variable treatgr as a factor variable and a regular variable is not supported

How can I change my code to avoid this error? I read the listers’ answers above, but do not know how to use them in my situation.


Thank you.

[Clyde Schechter]

The problem is in your -areg- command where treatgr appears both by itself with no prefix, and then in an interaction term where it is prefixed with i. When you specify a variable in a regression command with no prefix and outside of an interaction, Stata by default interprets it as c.-prefixed. But then you have it i.-prefixed right after that. Contradiction. Change that first treatgr to i.treatgr.

[Katherine Adams]

Clyde,

Thank you!

[Katherine Adams]

I simplified my model a little bit; in particular, now I use a variable i.treatalt instead of i.treatgroup:
areg lconsum i.treatalt tp i.treatalt#c.treat_numb_of_days, absorb(location) vce(cluster location)
margins, dydx(treatalt) at(treat_numb_of_days = (-403(86)285))
marginsplot

I got the following figure:


However, it is not the figure I want. I need to plot an event-study figure with energy consumption on the y-axis (I think it should be the interaction (‘i.treatalt#c.treat_numb_of_day’) coefficient, but I am not sure) and the number of days before/after the treatment on the x-axis (‘treat_numb_of_days’).

Does anyone know how can I plot my figure using “margins” command?

Thank you.

[Weiwen Ng]

I simplified my model a little bit; in particular, now I use a variable i.treatalt instead of i.treatgroup:
areg lconsum i.treatalt tp i.treatalt#c.treat_numb_of_days, absorb(location) vce(cluster location)
margins, dydx(treatalt) at(treat_numb_of_days = (-403(86)285))
marginsplot

I got the following figure:
[ATTACH=CONFIG]n1484732[/ATTACH]

However, it is not the figure I want. I need to plot an event-study figure with energy consumption on the y-axis (I think it should be the interaction (‘i.treatalt#c.treat_numb_of_day’) coefficient, but I am not sure) and the number of days before/after the treatment on the x-axis (‘treat_numb_of_days’).

Does anyone know how can I plot my figure using “margins” command?

Thank you.

Well, in your original code, it looks like you are plotting the (model-estimated) difference between treatment group 1 and the control group. That is one perfectly acceptable way to show the results.

If you wanted to plot the model-estimated energy consumption for each group, then I think your code should look something like:

[CODE]margins treatalt, at(treat_numb_of_days = (-403(86)285))[
marginsplot, xdimension(treat_numb_of_days)/CODE]

[Katherine Adams]

Weiwen,

Thank you! I will try using the code you specified.