rvaradhan commented 3 years ago

Suppose I have a 3-level categorical variable (groups labeled as A, B, C). I would like to compare the functional response (e.g., blood pressure) between A and B, and also A and C (i.e. A is my reference group). I would also like to adjust for potential confounders in the regression model.

Please see the attached code below that shows my attempt at this. I would like to know if this is a valid approach.

Thank you.

(0) Load dataset:

dataset = spm1d.data.uv1d.regress.SpeedGRF() Y,x = dataset.get_data() Ymax = np.max(Y, axis=1)

xp = x + 0.2*np.random.randn(x.size)

specify design matrix:

nCurves = x.size nFactors = 4 X = np.zeros((nCurves,nFactors)) X[:,0] = 1(Ymax <= 1.8) X[:,1] = 1(Ymax <= 2.1) & (Ymax > 1.8) X[:,2] = 1*(Ymax > 2.1) X[:,3] = xp # confounder to be adjusted for

specify a contrast vector:

c = [1,0,-1,0] # contrast between two-levels of a 3-level factor variable

(1) Conduct general linear test:

alpha = 0.05 roi = np.array(Y.shape[1][False]) roi[20:80] = [True]60

roi = np.array(roi, dtype=np.bool)

tstat = spm1d.stats.glm(Y, X, c, roi=roi) print(tstat.fwhm) # full width at half maximum = a measure of smoothness ti = tstat.inference(alpha, two_tailed=False, interp=True) print(ti.clusters)

(2) Plot:

plt.close('all') ax = plt.axes() ti.plot() ti.plot_threshold_label(fontsize=8) ti.plot_p_values(size=10, offsets=None) ax.set_xlabel('Time (%)') plt.show()

0todd0000 commented 3 years ago

The code looks fine in general. Whether or not this approach is valid depends on a variety of factors. For example:

Ymax appears to be a continuous variable, but is treated as a categorical variable in this example; this may or may not be valid.
The selected roi may or may not be justifiable in an a priori sense.
One-tailed inference may or may not be justifiable in an a priori sense.
If you only test the contrast [1,0,-1,0], then alpha is fine as is, but if you also wish to test other contrasts (e.g. [1,-1,0,0]) then a correction for multiple comparisons should be used. Although spm1d does not support these corrections directly, you could a Bonferroni correction as a conservative option, like this:

alpha      = 0.05 
ntests     = 2 
p_critical = spm1d.util.p_critical_bonf(alpha, ntests)                                     
ti         = tstat.inference(p_critical, two_tailed=False, interp=True)

rvaradhan commented 3 years ago

The construction of a categorical variable with 3-levels was purely artificial and only for illustrative purposes, hence the first comment is not a real issue. Also, the selection of ROI was also for illustrative purpose.

I appreciate the comment about Bonferroni correction for multiple contrasts.

Thank you.

From: Todd Pataky notifications@github.com Sent: Sunday, September 27, 2020 9:50:56 PM To: 0todd0000/spm1d Cc: Ravi Varadhan; Author Subject: Re: [0todd0000/spm1d] Effects of categorical variables in spm1d regression models (#145)

  External Email - Use Caution

The code looks fine in general. Whether or not this approach is valid depends on a variety of factors. For example:

Ymax appears to be a continuous variable, but is treated as a categorical variable in this example; this may or may not be valid.
The selected roi may or may not be justifiable in an a priori sense.
One-tailed inference may or may not be justifiable in an a priori sense.
If you only test the contrast [1,0,-1,0], then alpha is fine as is, but if you also wish to test other contrasts (e.g. [1,-1,0,0]) then a correction for multiple comparisons should be used. Although spm1d does not support these corrections directly, you could a Bonferroni correction as a conservative option, like this:

alpha = 0.05 ntests = 2 p_critical = spm1d.util.p_critical_bonf(alpha, ntests) ti = tstat.inference(p_critical, two_tailed=False, interp=True)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2F0todd0000%2Fspm1d%2Fissues%2F145%23issuecomment-699726709&data=02%7C01%7Cravi.varadhan%40jhu.edu%7Cf27e6a86efab44310e0608d86350f1ed%7C9fa4f438b1e6473b803f86f8aedf0dec%7C0%7C0%7C637368546619675898&sdata=7cdp3KhQRyvklENWPw3JLCpcKjySlpribj7MRZWN8Oo%3D&reserved=0, or unsubscribehttps://nam02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAJZCQA2G7FAXCYDXBIJHVLSH7TYBANCNFSM4R3VAPHQ&data=02%7C01%7Cravi.varadhan%40jhu.edu%7Cf27e6a86efab44310e0608d86350f1ed%7C9fa4f438b1e6473b803f86f8aedf0dec%7C0%7C0%7C637368546619675898&sdata=uJJEH2zp7IAwpXuD%2FHmIQlrbw0q8yTy167oHWZqiY2w%3D&reserved=0.