[William Lisowski]

Added in edit: Post #15 was updated while I was writing this post. But I'm still out of here until tomorrow my time.

Refresh your memory of post #7 above.

Then, please reflect on the advice of the Statalist FAQ linked to from the top of the page, as well as from the Advice on Posting link on the page you used to create your post.

12. What should I say about the commands and data I use?

Help us to help you by producing self-contained questions with reproducible examples that explain your data, your code, and your problem. This helps yet others too, as they will find it easier to learn from your questions and the answers to them.

Perhaps if you accompany your example data with the entire code you used - the corrected version of the swrct program and whatever other code was in your do-file - so that someone can start by reading in your data from post #15 and then run exactly the code you show to reproduce the error, someone will be able to help.

Or perhaps someone will guess how variables in your data from post #15 relate to the variables in your code in post #9, and will be interested in modifying the code in post #9 to incorporate the correction I pointed out and any other changes needed to get it to run on the data from post #15.

I'm out of here for the night, so you may again have lost yourself a chance at a quick answer. But Joro's profile tells us he's in London, so while we haven't seen him in the last few hours, perhaps with adequate information he will be interested in trying again while I'm away from Statalist.

[CEdward]

Thank you William.

So this is the code that I am running which gives that error r(509); described in #15. What I think is perhaps the problem is something to do with the xtgeebcv command itself as this also happens when I run the sample data and code provided in the installation.

I hope I was clear in articulating the relevant details.

clear *
set more on

*Generate multi-level data
capture program drop swcrt
program define swcrt, rclass
args num_clus clussize intercept intrvcoeff timecoeff1 timecoeff2 timecoeff3 timecoeff4 timecoeff5 sigma_u3 sigma_error alpha
//display "num_clus = `num_clus' and clussize = `clussize' and intercept = `intercept' and timecoeff1 = `timecoeff1' and timecoeff2 = `timecoeff2' and timecoeff3 = `timecoeff3' and timecoeff4 = `timecoeff4' and timecoeff5 = `timecoeff5' and sigma_u3 = `sigma_u3' and sigma_error = `sigma_error' and alpha = `alpha'"
assert `num_clus' > 0 & `clussize' > 0 & `intercept' > 0 & `intrvcoeff' > 0 & `timecoeff1' < 0 & `timecoeff2' < 0 & `timecoeff3' < 0 & `timecoeff4' < 0 & `timecoeff5' < 0 & `sigma_u3' > 0 & `sigma_error' > 0 & `alpha' > 0
*Generate simulated multi—level data
qui
clear
set obs `num_clus'
qui gen cluster = _n
qui gen group = 1+mod(_n-1,4)
*Generate cluster-level errors
qui gen u_3 = rnormal(0,`sigma_u3')
expand `clussize'
bysort cluster: gen individual = _n
*Set up time
expand 6
bysort cluster individual: gen time = _n-1
*Set up intervention variable
gen intrv = (time>=group)
*Generate residual errors
qui gen error = rnormal(0,`sigma_error')
*Generate outcome y
qui gen y = `intercept' + `intrvcoeff'*intrv + `timecoeff1'*1.time + `timecoeff2'*2.time + `timecoeff3'*3.time + `timecoeff4'*4.time + `timecoeff5'*5.time + u_3 + error

*Fit multi-level model to simulated dataset
*xtset cluster
*xtgee y intrv i.time, family(gaussian) link(identity) corr(exchangeable)

xtgeebcv y intrv i.time, cluster(cluster) family(gaussian) link(identity) corr(exchangeable) stderr(kc)
*Return estimated effect size, bias, p-value, and significance dichotomy
tempname M
matrix `M' = r(table)
return scalar b = _b[intrv]
return scalar bias = _b[intrv] - `intrvcoeff'
return scalar p = `M'[1,4]
return scalar p_= (`M'[1,4] < `alpha')
exit
end swcrt

*Stepped-wedge specifications
local num_clus 3 6 9 18 36
local clussize 5 10 15 20 25
*Model specifications
local intercept 17.87
local timecoeff1 -5.42
local timecoeff2 -5.72
local timecoeff3 -7.03
local timecoeff4 -6.13
local timecoeff5 -9.13
local intrvcoeff 5.00
local sigma_u3 25.77
local sigma_u2 120.62
local sigma_error 38.35
local alpha 0.05
local nrep 1

*Postfile to store results
tempname step
tempfile powerresults
capture postutil clear
postfile `step' num_clus clussize intrvcoeff b p p_ bias using `powerresults.dta', replace
*Loop over number of clusters
foreach c of local num_clus{
display as text "Number of clusters" as result "`c'"
foreach s of local clussize {
display as text "Cluster size" as result "`s'"
forvalue i = 1/`nrep'{
display as text "Iterations" as result `nrep'
quietly swcrt `c' `s' `intercept' `intrvcoeff' `timecoeff1' `timecoeff2' `timecoeff3' `timecoeff4' `timecoeff5' `sigma_u3' `sigma_error' `alpha'
post `step' (`c') (`s') (`intrvcoeff') (`r(b)') (`r(p)') (`r(p_)') (`r(bias)')
}

}

}
postclose `step'

*Open results, calculate power
use `powerresults', clear

levelsof num_clus, local(num_clus)
levelsof clussize, local(clussize)
matrix drop _all
*Loop over combinations of clusters
*Add power results to matrix
foreach c of local num_clus{
foreach s of local clussize{
quietly ci proportions p_ if num_clus == `c' & clussize = `s'
local power `r(proportion)'
local power_lb `r(lb)'
local power_ub `r(ub)'
quietly ci mean bias if num_clus == `c' & clussize = `s'
local bias `r(mean)'
matrix M = nullmat(M) \ (`c', `s', `intrvcoeff', `power', `power_lb', `power_ub', `bias')
}
}

*Display the matrix
matrix colnames M = c s intrvcoeff power power_lb power_ub bias //Arguments should match those in the previous line to avoid conformability errors
matrix list M, noheader format(%3.2f)

[Joro Kolev]

I can confirm that there is a bug in the user written -xtgeebcv- , it produces an error on authors example dataset. This is the statement that generates the error:

= local name: word 3 of "Robust" "Degrees-of-Freedom" "Kauermann-Carroll" "Mancl-DeRouen" "Fay-Graub
> ard" "Morel-Bokossa-Neerchal"
- if "`stderr'" == "`bcv'" {
= if "kc" == "kc" {
- local dof=nclust[1,1]-ncoef[1,1]
- local last = rowsof(e(V))+1
matrix operators that return matrices not allowed in this context

. use mkvtrial.dta, clear

. xtgeebcv know i.arm i.stratum i.ethnicgp, family(binomial) link(logit) cluster(community) stderr(k
> c) nolog
 
Note: Family is binomial and link is logit
Using exchangeable working correlation
with scale parameter divided by K - p
       panel variable:  community (unbalanced)

GEE population-averaged model                   Number of obs     =      4,100
Group variable:                  community      Number of groups  =         20
Link:                                logit      Obs per group:
Family:                           binomial                    min =        169
Correlation:                  exchangeable                    avg =      205.0
                                                              max =        257
                                                Wald chi2(4)      =      43.75
Scale parameter:                         1      Prob > chi2       =     0.0000

------------------------------------------------------------------------------
        know |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       1.arm |   .8270697   .1482322     5.58   0.000     .5365398    1.117599
             |
     stratum |
          2  |   .0503953    .179305     0.28   0.779    -.3010361    .4018266
          3  |   .1252695   .1924411     0.65   0.515    -.2519082    .5024471
             |
  1.ethnicgp |  -.3040093   .0879244    -3.46   0.001    -.4763379   -.1316807
       _cons |  -.0108448    .164238    -0.07   0.947    -.3327454    .3110558
------------------------------------------------------------------------------
matrix operators that return matrices not allowed in this context
r(509);

.

[Joro Kolev]

What happens is what I explained in #10: Stata does not allow matrix expressions that require evaluation to be equated to scalars/locals, even if those matrix expressions evaluate to scalars.

Here is the problem with the authors' code replicated on the auto data:

. sysuse auto, clear
(1978 Automobile Data)

. qui reg price mpg headroom

. local rowsplusone = rowsof(e(V))+1
matrix operators that return matrices not allowed in this context
r(509);

[Joro Kolev]

The fix to this problem is to split the operation in two rows, two steps. This is very annoying, but I do not know of any other solution. Here:

. mat tempmat =  rowsof(e(V))+1

. local rowsplusone = tempmat[1,1]

What happens here is that the first operation which requires matrix evaluation, we set equal to a matrix (although the result of the matrix evaluation is equal to a scalar). Then in the second step we simply pick up the matrix element and equate it to the local, and for this latter Stata does not complain.

[Joro Kolev]

Bjarte Aagnes has resolved the mystery of how the authors of -xtgeebcv- can possibly not notice that their command is generating errors on the example dataset.

The problem described in #18, arises in Stata 15.1, and I presume in all previous versions of Stata, but does not arise in Stata 16.1.

[William Lisowski]

Following on from the discoveries from Joro Kolev and Bjarte Aagnes, there are two available versions of xtgeebcv

• Version 2.1 24Jan2020 is distributed from Stata Journal as package sj0599
• Version 3.0 04Mar2020 is distributed from SSC as package xtgeebcv

The problem code appears only in version 3.0.

Both versions include version control for version 15. However, I'd guess that version 3.0 was only tested on version 16, without version control to the earlier version, as Bjarte suggests.

Perhaps John Gallis or Fan Li will see this and update the SSC version.

Added in edit: version control is more subtle than my casual knowledge led me to understand.The following example was run on an installation of Stata 16.1, and is a reproducible example of the problem encountered with xtgeebcv. Despite the appearance of version 15 control in the program definition, it runs correctly on version 16.1 unless version control is additionally enforced when the command is called.

. sysuse auto, clear
(1978 Automobile Data)

. program define gnxl
  1. version 15
  2. quietly regress price mpg headroom
  3. local rowsplusone = rowsof(e(V))+1
  4. display "rowsplusone is `rowsplusone'"
  5. end

. gnxl
rowsplusone is 4

. version 15: gnxl
matrix operators that return matrices not allowed in this context
r(509);

[CEdward]

The fix to this problem is to split the operation in two rows, two steps. This is very annoying, but I do not know of any other solution. Here:

. mat tempmat = rowsof(e(V))+1

. local rowsplusone = tempmat[1,1]

What happens here is that the first operation which requires matrix evaluation, we set equal to a matrix (although the result of the matrix evaluation is equal to a scalar). Then in the second step we simply pick up the matrix element and equate it to the local, and for this latter Stata does not complain.

Joro, thank you for your input and investigation. William, yourself as well.

I had suspected that it is a version control issue as I noticed that previously that William was able to run it just fine. Indeed, I am using version 15.1.

I follow your fix to the problem.

So it seems to me that unless the author's fix the bug then this is not usable for me because my program will break each time. In this instance, is it common to request for the authors to update?

[William Lisowski]

The which and doedit commands demonstrated below allow you to

• open the installed version of xtgeebcv.ado in your do-file editor
• where you can make the changes described in post #20
  • I would comment out the faulty command and add the two recommended commands after it, just in case you need to undo your changes
• then save the corrected xtgeebcv.ado
• then from the Stata Command window run clear all to ensure that the next time you run xtgeebcv the (corrected) ado file will be loaded

.  which xtgeebcv
/Users/lisowskiw/Library/Application Support/Stata/ado/plus/x/xtgeebcv.ado
*!version3.0 04Mar2020

. doedit "/Users/lisowskiw/Library/Application Support/Stata/ado/plus/x/xtgeebcv.ado"

Added in edit: a one-line solution

doedit "`c(sysdir_plus)'/x/xtgeebcv.ado"

[CEdward]

Thank you to William and Joro for their help on solving this problem.

For other Stata users, this thread had two problems: incorrect calling of the program and a bug in the xtgeebcv package that was not updated by the developers.

The solution to the latter is found in posts #20 and #24.

[John Gallis]

Hi all. In the future, please email me with bugs in my programs, as I am not alerted by email with posts on statalist, and I rarely check here. Please note that the version of xtgeebcv on Stata Journal has been tested and should be bug-free. I added some functionality to the SSC version, and I will look into this bug. It appears that something I added in the 3.0 version only works on version 16. I will update the code on SSC to only allow running the program on version 16. Thanks!

[William Lisowski]

Post #26 leads me to mention that Statalist members can configure the Notifications tab in their Profile to (a) send email to them when they are mentioned in a Statalist post, as in post #22 above, and (b) when they are sent a Private Message on Statalist, while (c) allowing them to turn off all the other options for email. Statalist members who have provided community contributed software have found that a helpful way of engaging with their user community through Statalist at a minimal increase email.

[John Gallis]

Thanks, William Lisowski! I will set that up now, as that will be very helpful (although I have found sometimes in the past my programs are being discussed on Statalist and no one tags me, so thanks for tagging me in this one).

Also, that was fast, but now an updated version of xtgeebcv has been uploaded to SSC which only allows running on version 16.1 or later, and the help file directs those with version 15 to download the Stata Journal program. Once Stata 17 is released, I will probably update the Stata Journal version to match the SSC version to avoid confusion. Hopefully by then most people still using version 15 will at least have updated to version 16.