problem with bcoaff

[Carolina Gracitelli]

Hi,

Sorry to disturb. I'm a fellow of Department of Ophthalmology, UCSD, USA. Recently I installed the command of bcoeff written by you many years ago to calculate the loss rate of retinal nerve fiber layer thickness over time in glaucoma patients. However, I found something weird: why did I get so differential results with similar raw data? I'm not so sure the calculating details of bcoeff, so if you don't mind, could you please take a look at my dataset and tell me the reason? I do appreciate your kind help.

Regards,

Carolina






(For example, SD1031 had both eyes with the same follow-up time and visitnumber, but the bcoeffrnfl of his right eye and left eye was so different. The same question exited for SD2011, SD2174, SD2203��
My codes were: "bcoeff avgthickness time, by(id) g(bcoeffrnfl)", in which avgthickness was the same to rnfl_g except for variable type and id was generated with maskedid )

[Nick Cox]

"You" is not clear here, but I looked up bcoeff and found that it is on SSC and written by Zhiqiang Wang and myself. The program is dated 20 June 2000 and the 2004 advice to use other stuff is already a decade old. (Last time I looked Zhiqiang had left ANU, so that email address will bounce.) Regardless of that, I see no way to comment on your dataset.

----------------------------------------------------------------------------------------
package bcoeff from http://fmwww.bc.edu/repec/bocode/b
----------------------------------------------------------------------------------------

TITLE
'BCOEFF': module to save regression coefficients to new variable

DESCRIPTION/AUTHOR(S)

bcoeff saves in a new variable regression coefficients (more
generally, the b coefficient from a regression-like model) for
each of several groups of observations. (bcoeff supersedes
deltaco by Zhiqiang Wang.) Note added 2004: this package was
written for Stata 6 and is now superseded for most purposes by
use of foreach or statsby in Stata 7 or 8.

Author: Zhiqiang Wang, Menzies School of Health Research
Support: email [email protected]

Author: Nicholas J. Cox, University of Durham
Support: email [email protected]

Distribution-Date: 20000620

[Ting Liu]

Thank you, Nick. The question is actually mine, but I didn't have the access to post it yesterday. Now my registration was verified, so I'm available to take it over. This is the screen shots of my dataset, and hope they can be clear and helpful.

Thanks again.

[Nick Cox]

I don't take anything from this yet except that you find your results puzzling. As you are not showing any results, I don't know what sort of comment you think is possible.

[Ting Liu]

As you can see from the highlighted area, SD1031 had both eyes with the same follow-up time and visitnumber, and similar "rnfl_g" data over time, but the "bcoeffrnfl" of his right eye and left eye was so different. The same question exited for SD2011, SD2174...
My codes were: "bcoeff avgthickness time, by(id) g(bcoeffrnfl)", in which avgthickness was the same value to rnfl_g except for variable type and id was generated with maskedid.

Thank you, Nick.

[Nick Cox]

That program is just a reporting program. Nothing you've reported so far implies that it is mangling results. If you suspect that, you should run the regressions separately and check that results are the same (or show the difference if there is one). You are slowly adding details but I'm afraid that the flavour remains that you are puzzled by the results. You have a clear course of action: plotting the two or more regressions to see why they differ. That's a matter of scatter plots and regression lines. I don't have access to Stata right now but in any case screenshots of the data don't allow easy checking of anything.

[Ting Liu]

Sorry for slowing adding the details and thanks again for the patience with a new comer. However, I didn't find the difference between the right eye's rnfl_g and left eye's rnfl_g with the same maskedid in the scatterplot so big as the ones in bcoeffrnfl. Frankly, I'm indeed confused.
Since I can't upload the dataset as .dta or .csv file, if you could spare me a few minutes, the raw dataset can be downloaded from my Dropbox: https://www.dropbox.com/home/For%20Nick. I do appreciate your help and take your time. Have a nice weekend.

[Joseph Coveney]

Why can't you attach the dataset to your post?

If I've copied the data below for Patient SD1031 correctly from your screenshot attachments, then I don't understand why you have trouble with the finding that the regression coefficients are so different between eyes for this patient.

version 13.1

clear *
set more off
input double year int thickness_r thickness_l
0        101  98
1.034908  99 101
1.533916  98  99
2.069815  98  99
2.792608  96 101
3.08282   97  99
end

sort year
list, noobs separator(0) abbreviate(20)

regress thickness_r c.year
regress thickness_l c.year

graph twoway lfit thickness_r year, lcolor(blue) || ///
    lfit thickness_l year, lcolor(red) || ///
    scatter thickness_r year, mcolor(blue) || ///
    scatter thickness_l year, mcolor(red) ///
        ylabel(95(2.5)102.5, angle(horizontal) nogrid) ///
        ytitle(Retinal Nerve-fiber Layer Thickness) ///
        xtitle(Time (y)) legend(off)

exit

From the values in your screenshot (again, assuming that I've correctly transcribed them), I get vastly different slopes than you get using avgthickness.

[Ting Liu]

Thanks a lot, Joseph. Yes, the slopes from your codes are vastly different from mine, which was got using command "bcoeff".
Could you tell me how I can upload the dataset? It seems unaccessible with the function of "upload attachment", which always tells me: "wrong format..."

[Joseph Coveney]

I suspect that the difference between the slopes that you and I got has nothing to do with bcoeff and everything to do with avgthickness. Regardless, I still don't see why you have a problem with the left-right difference for Patient SD1031 (and presumably for the others, as well).

I have the same problem attaching a Stata dataset (.dta) file as you do. I get the error message, "This is not a valid image file. SD1031.dta", which should probably be brought to the attention of the administrators. But you can save the data as a text file in Stata's dictionary format using outfile using <filename>.txt, dictionary and that can be attached. At least I could with Patient SD1031's transcribed data (see attachment below).

[Ting Liu]

Thank you, Joseph. But the "rnfl_g" has the same value to the "avgthickness" except for variable format. Now I rename the variable name of "rnfl_g" into "avgthickness", please check my Dropbox for the raw dataset "https://www.dropbox.com/home/For%20Joseph". Sorry, the dataset exceeds the up limit of attachment.
Thanks again.

[Joseph Coveney]

As you can see from the highlighted area, SD1031 had both eyes with the same follow-up time and visitnumber, and similar "rnfl_g" data over time, but the "bcoeffrnfl" of his right eye and left eye was so different. The same question exited for SD2011, SD2174...
My codes were: "bcoeff avgthickness time, by(id) g(bcoeffrnfl)", in which avgthickness was the same value to rnfl_g except for variable type and id was generated with maskedid.

Ting, think about it for a moment: according to the highlighted area for Patient SD1031 oculus dexter, the average thickness gradually declines from 101 units to 97 units (a change of -4 units of thickness) over a period of a little over three years, yet your dataset shows a slope coefficient of +21.98576 (units per year). I don't know what your complete code actually was, but it couldn't possibly have been what you claim.

Here, I downloaded your dataset from dropbox.com, cleaned up the UTF-8 Unicode that it contained, renamed avgthickness back to rnfl_g, and ran the following do-file. You'll see that there is nothing strange about the user-written bcoeff compared to what I got with regress manually above.

.version13.1

.
.clear*

.setmoreoff

.
.useTing

.
.renamebcoeffrnflbcoeffrnfl_old

.destringrnfl_g,generate(avgthickness)
rnfl_ghasallcharactersnumeric;avgthicknessgeneratedasint

.
.egenintid=group(maskedideye)

.bcoeffavgthicknesstime,by(id)g(bcoeffrnfl)

.
.sortmaskedideyevisitnumber

.listmaskedideyebcoe*if///
>inlist(maskedid,"SD1031":maskedid00,"SD2011":maskedid00,///
>"SD2174":maskedid00)&!time,noobssepby(maskedid)abbreviate(20)

+----------------------------------------------+
|maskedideyebcoeffrnfl_oldbcoeffrnfl|
|----------------------------------------------|
|SD1031OD21.98576-1.424926|
|SD1031OS-3.952217.3641296|
|----------------------------------------------|
|SD2011OD-1.696892-1.696892|
|SD2011OS8.184894-.8293823|
|----------------------------------------------|
|SD2174OD13.24641-1.233431|
|SD2174OS-1.298882-1.298882|
+----------------------------------------------+

.
./*encodeeye,generate(eye_c)label(Eye)
>mixedavgc.time##i.eye_c||maskedid:time,covariance(unstructured)||///
>eye_c:time,covariance(unstructured)remlstddeviationsnolrtestnolog
>
>predictdoubleresiduals,residuals
>predictdoublefitted,fitted
>graphtwowayscatterresidualsfitted,mcolor(black)msize(small)///
>yline(0,lcolor(gs8))ylabel(,angle(horizontal)nogrid)
>
>qnormresiduals,ylabel(,angle(horizontal)nogrid)
>pnormresiduals,ylabel(,angle(horizontal)nogrid)*/
.
.exit

endofdo-file


.

What you got before (bcoeffrnfl_old) is only occasionally the same as what I get with bcoeff and often very different on your own dataset.

I've zipped the dataset and my do-file, renamed the zip file's file-name suffix as .png in order to be able to attach it (contrary to what the forum's Software FAQ says, zip files are not allowed), and attached it below. You can run the do-file for yourself on your dataset, by saving the attachment, renaming it back and unzipping it on your local disc.

[Nick Cox]

As bcoeff is partly my program, and its first author is not in sight, I retain an interest in this question. Joseph picked it up, I guess out of curiosity mixed with altruism, while I was travelling. The exchanges since I last commented haven't changed my impression that whatever problem there is here is nothing to do with the use of bcoeff. As soon as Ting posts some results that clearly imply otherwise, I will try a closer look.

Meanwhile I'll repeat a comment made a while back. bcoeff remains accessible so that people can replicate anything done with it but already in 2004 people were advised to use something else. That's not because I have reason to doubt bcoeff, just that changes in Stata have made it redundant.

[Ting Liu]

1) Thank you both so much for the altruism and the great help in such kind of detail. Yes, There're only 2 patients (SD4022 and SD4107) with the same results for bcoeffrnfl_old and bcoeffrnfl and I've double checked the raw data on the machine to ensure the validity and reasonability.

Notice there're some huge difference between the bcoeffrnfl_old and bcoeffrnfl in the same patient and the same eye. Probably it's because what Nick told us that "changes in Stata have made it redundant".

2) Now I'm a little concerned about another question: will the command xtmixed produces some results like bcoeff? Actually, I just used bcoeff to double check the outliers of the results from xtmixed and happened to get the above weird findings. However, I finally must return to use the linear mixed model to calculate the rate of change of avgthickness, considering other confounding factors, such as age, race… This is a part of my do file to calculate the rate of change but did not include the confounding factors yet:

use "/Users/air/Documents/Stata/Intereye/dta/dropped some outliers and counted NASDF as progressors_IOP_BY.dta", clear

bysort maskedid eye: gen id=1 if visitnumber==1
replace id=sum(id)

xtmixed avgthickness time if PROG_OD==1 | PROG_OD==2 || id: time, cov(un)
predict predictedrnfl, fitted
xtline predictedrnfl if PROG_OD==1 | PROG_OD==2, t(time) i(id) overlay

xtmixed avgthickness time if PROG_OD==3 || id: time, cov(un)
predict predictedrnfl, fitted
xtline predictedrnfl if PROG_OD==3, t(time) i(id) overlay

xtmixed avgthickness time if PROG_OD==4 | PROG_OD==5 || id: time, cov(un)
predict predictedrnfl, fitted
xtline predictedrnfl if PROG_OD==4 | PROG_OD==5, t(time) i(id) overlay

/* PROG_OD==1: progrssed in right eye; PROG_OD==2: progrssed in left eye; PROG_OD==3: fellow eye of the progressor; PROG_OD==4: right eye of the nonprogressor; PROG_OD==5: left eye of the nonprogressor */

As was shown in the previous work by another fellow, it is supposed to be of significant loss of avgthickness in PROG_OD==1/2/3, and even in PROG_OD==4/5, who is the nonprogressor, there supposes to be of physiological loss of avgthickness. However, I didn't get the expected results either. Is there something wrong with my dataset or the harmony between STATA and the command?
Many thanks in advance.

[Ting Liu]

The attached are some results and graphs of the above do file for xtmixed:






If you can spare some time to check the dataset for the final analysis, please refer to "https://www.dropbox.com/home/For%20Nick%20and%20Joseph?d=1".
Thanks again.