Untitled

Just ran a moderation analysis using the NESARC dataset to see if marital status changes the relationship between personal income and number of children under 18 in the household.

Here’s the TL;DR:

💡 Main finding: Income is (still) barely related to number of kids. BUT—if you’re married, the pattern shifts just a little.

📉 For non-married folks: Higher income is very slightly linked to fewer kids. (We’re talking fractions of a child per $10,000...)

📈 For married folks: That negative trend basically flattens out. So among married people, income doesn’t really predict number of children at all.

👫 Also: Married people, regardless of income, tend to have more kids than non-married folks. No surprise there, but it’s nice to see the data back it up.

📊 Here's the math-y part (if you're into that): We ran a regression model including:

Income (S1Q10A)

Marital status (MARITAL)

And the interaction term (Income × Married)

📍The interaction was statistically significant (p < .001), but the effect size? Tiny. Model only explained 4.2% of the variance in number of kids. (So yeah, life is more complicated than one variable!)

🤷‍♀️ In short: Marriage nudges the income-kid link around a bit, but money still isn’t a strong predictor of how many kids you have—married or not.

So if you’re thinking: "Do rich people have fewer kids because they’re rich?" or "Does marriage make income matter more for family size?"

The data kind of shrugs and says: "Not really."

In case you want to try out the code: ----------------------------

import pandas as pd import numpy as np import seaborn as sns import scipy.stats import matplotlib.pyplot as plt import statsmodels.formula.api as smf

Load the data

data = pd.read_csv('nesarc.csv', low_memory=False)

Convert to numeric properly

data['CHLD0_17'] = pd.to_numeric(data['CHLD0_17'], errors='coerce') data['S1Q10A'] = pd.to_numeric(data['S1Q10A'], errors='coerce') data['MARITAL'] = pd.to_numeric(data['MARITAL'], errors='coerce')

Drop rows with NA in either of the two columns

data_clean = data[['S1Q10A', 'CHLD0_17', 'MARITAL']].dropna()

data_clean['Married'] = data_clean['MARITAL'].apply(lambda x: 1 if x == 1 else 0)

Scatterplot

sns.lmplot(x='S1Q10A', y='CHLD0_17', hue='Married', data=data_clean, palette='Set1', height=6, aspect=1.2) plt.xlabel('Total Personal Income (USD)') plt.ylabel('Number of Children Under 18') plt.title('Income vs Number of Children by Marital Status') plt.show()

Fit linear regression with interaction term

model = smf.ols('CHLD0_17 ~ S1Q10A * Married', data=data_clean).fit()

Output regression summary

print(model.summary())

plt.figure(figsize=(8,6)) sns.regplot(data=data_clean[data_clean['Married'] == 1], x='S1Q10A', y='CHLD0_17', label='Married', scatter_kws={'alpha':0.3}) sns.regplot(data=data_clean[data_clean['Married'] == 0], x='S1Q10A', y='CHLD0_17', label='Not Married', scatter_kws={'alpha':0.3}, color='green') plt.xlabel('Total Personal Income (USD)') plt.ylabel('Number of Children Under 18') plt.title('Income vs Number of Children by Marital Status') plt.legend() plt.tight_layout() plt.show()

Just ran a Pearson correlation between total personal income (S1Q10A) and number of children under 18 in the household (CHLD0_17) using the NESARC dataset. The result?

r = -0.016 p = 0.0009

💥 Yes, it's statistically significant. 😅 But practically? Not so much.

This teeny-tiny negative correlation means higher income is slightly linked to fewer kids—but the effect is so small it’s basically a shrug from the universe. 🤷‍♂️📉

In short: Money ≠ more (or fewer) children — at least not in any meaningful way.

in case you want to try out the code:

import pandas as pd import numpy as np import seaborn as sns import scipy.stats import matplotlib.pyplot as plt

Load the data

data = pd.read_csv('nesarc.csv', low_memory=False)

Convert to numeric properly

data['CHLD0_17'] = pd.to_numeric(data['CHLD0_17'], errors='coerce') data['S1Q10A'] = pd.to_numeric(data['S1Q10A'], errors='coerce')

Drop rows with NA in either of the two columns

data_clean = data[['S1Q10A', 'CHLD0_17']].dropna()

Scatterplot

sns.regplot(x="S1Q10A", y="CHLD0_17", fit_reg=True, data=data_clean) plt.xlabel('Total Personal Income (USD)') plt.ylabel('Number of Children Under 18 in Household') plt.title('Scatterplot: Income vs. Number of Children') plt.show()

Pearson correlation

Using data from the NESARC survey, we dug into how your relationship status might influence how much you drink. The stats? Unfiltered and telling.

📊 Key Insight:

There’s a significant association between marital status and alcohol consumption (Chi² = 392.18, p < .001).

Married individuals top the charts across all drinking levels.

Divorced/Separated drink moderately — but steadily.

Never Married lean toward lighter drinking, but mid-range isn't rare.

Widowed appear quietly in higher consumption brackets. (Coping, perhaps?)

⚖️ Pairwise Comparisons:

We broke it down further with post-hoc testing:

Married vs Never Married → Big difference (p < .001)

Married vs Divorced/Cohabiting → Still distinct (p < .01)

Never Married vs Separated → Subtle, yet significant (p ≈ .015)

Your marital status might be quietly shaping your drinking habits — whether you're toasting to love, numbing the breakup, or just navigating life solo. In case you want to try the code yourself ````````````````````````

import pandas as pd import numpy as np import scipy.stats as stats import seaborn as sns import matplotlib.pyplot as plt

Load the dataset

data = pd.read_csv('nesarc.csv', low_memory=False)

Convert variables to numeric

data['S2AQ8A'] = pd.to_numeric(data['S2AQ8A'], errors='coerce') # Alcohol frequency data['MARITAL'] = pd.to_numeric(data['MARITAL'], errors='coerce') # Marital status

Subset data (assuming sub1 already exists — else, create one)

sub1 = data.copy() sub2 = sub1.copy()

Replace invalid or missing codes with NaN

sub2['S2AQ8A'].replace([99, np.nan], np.nan, inplace=True) sub2['MARITAL'].replace(9, np.nan, inplace=True)

Filter only valid values

sub2 = sub2[(sub2['S2AQ8A'].between(1, 10)) & (sub2['MARITAL'].between(1, 6))]

Create a label-friendly MARITALSTATUS column (you can use actual labels instead of scores)

marital_map = { 1: 'Married', 2: 'Cohabiting', 3: 'Widowed', 4: 'Divorced', 5: 'Separated', 6: 'Never Married' } sub2['MARITALSTATUS'] = sub2['MARITAL'].map(marital_map)

-------------------------

🔍 Main Chi-Square Test

-------------------------

ct_main = pd.crosstab(sub2['S2AQ8A'], sub2['MARITALSTATUS']) print("Main Contingency Table:\n", ct_main)

print("\nColumn Percentages:\n", ct_main / ct_main.sum(axis=0))

chi2, p, dof, expected = stats.chi2_contingency(ct_main) print("\nChi-square test results:") print(f"Chi2 Value: {chi2:.2f}") print(f"p-value: {p:.4f}") print(f"Degrees of Freedom: {dof}") print("Expected Counts:\n", pd.DataFrame(expected, index=ct_main.index, columns=ct_main.columns))

-------------------------

📊 Visualization

-------------------------

plt.figure(figsize=(12, 6)) sns.countplot(data=sub2, x='MARITALSTATUS', hue='S2AQ8A', palette='Set2') plt.title("Drinking Frequency by Marital Status") plt.xlabel("Marital Status") plt.ylabel("Count of Respondents") plt.legend(title="Drinking Frequency Code", bbox_to_anchor=(1.05, 1), loc='upper left') plt.tight_layout() plt.show()

-------------------------

🔁 Optional: Pairwise Comparisons

-------------------------

def chi_square_pairwise(df, group_col, compare_col, groups): sub_df = df[df[group_col].isin(groups)] ct = pd.crosstab(sub_df[compare_col], sub_df[group_col]) chi2, p, dof, expected = stats.chi2_contingency(ct) print(f"\nChi-square for {groups[0]} vs {groups[1]}") print(f"Chi2: {chi2:.2f}, p-value: {p:.4f}, DOF: {dof}") print("Contingency Table:\n", ct)

Run pairwise tests

pairwise_groups = ('Married', 'Never Married'), ('Married', 'Divorced'), ('Married', 'Cohabiting'), ('Never Married', 'Separated')

for g1, g2 in pairwise_groups: chi_square_pairwise(sub2, 'MARITALSTATUS', 'S2AQ8A', [g1, g2])

Using data from the NESARC survey, we dug into how your relationship status might influence how much you drink. The stats? Unfiltered and telling.

📊 Key Insight:

There’s a significant association between marital status and alcohol consumption (Chi² = 392.18, p < .001).

Married individuals top the charts across all drinking levels.

Divorced/Separated drink moderately — but steadily.

Never Married lean toward lighter drinking, but mid-range isn't rare.

Widowed appear quietly in higher consumption brackets. (Coping, perhaps?)

⚖️ Pairwise Comparisons:

We broke it down further with post-hoc testing:

Married vs Never Married → Big difference (p < .001)

Married vs Divorced/Cohabiting → Still distinct (p < .01)

Never Married vs Separated → Subtle, yet significant (p ≈ .015)

Your marital status might be quietly shaping your drinking habits — whether you're toasting to love, numbing the breakup, or just navigating life solo. In case you want to try the code yourself ````````````````````````

import pandas as pd import numpy as np import scipy.stats as stats import seaborn as sns import matplotlib.pyplot as plt

Load the dataset

data = pd.read_csv('nesarc.csv', low_memory=False)

Convert variables to numeric

data['S2AQ8A'] = pd.to_numeric(data['S2AQ8A'], errors='coerce') # Alcohol frequency data['MARITAL'] = pd.to_numeric(data['MARITAL'], errors='coerce') # Marital status

Subset data (assuming sub1 already exists — else, create one)

sub1 = data.copy() sub2 = sub1.copy()

Replace invalid or missing codes with NaN

sub2['S2AQ8A'].replace([99, np.nan], np.nan, inplace=True) sub2['MARITAL'].replace(9, np.nan, inplace=True)

Filter only valid values

sub2 = sub2[(sub2['S2AQ8A'].between(1, 10)) & (sub2['MARITAL'].between(1, 6))]

Create a label-friendly MARITALSTATUS column (you can use actual labels instead of scores)

marital_map = { 1: 'Married', 2: 'Cohabiting', 3: 'Widowed', 4: 'Divorced', 5: 'Separated', 6: 'Never Married' } sub2['MARITALSTATUS'] = sub2['MARITAL'].map(marital_map)

-------------------------

🔍 Main Chi-Square Test

-------------------------

ct_main = pd.crosstab(sub2['S2AQ8A'], sub2['MARITALSTATUS']) print("Main Contingency Table:\n", ct_main)

print("\nColumn Percentages:\n", ct_main / ct_main.sum(axis=0))

chi2, p, dof, expected = stats.chi2_contingency(ct_main) print("\nChi-square test results:") print(f"Chi2 Value: {chi2:.2f}") print(f"p-value: {p:.4f}") print(f"Degrees of Freedom: {dof}") print("Expected Counts:\n", pd.DataFrame(expected, index=ct_main.index, columns=ct_main.columns))

-------------------------

📊 Visualization

-------------------------

plt.figure(figsize=(12, 6)) sns.countplot(data=sub2, x='MARITALSTATUS', hue='S2AQ8A', palette='Set2') plt.title("Drinking Frequency by Marital Status") plt.xlabel("Marital Status") plt.ylabel("Count of Respondents") plt.legend(title="Drinking Frequency Code", bbox_to_anchor=(1.05, 1), loc='upper left') plt.tight_layout() plt.show()

-------------------------

🔁 Optional: Pairwise Comparisons

-------------------------

def chi_square_pairwise(df, group_col, compare_col, groups): sub_df = df[df[group_col].isin(groups)] ct = pd.crosstab(sub_df[compare_col], sub_df[group_col]) chi2, p, dof, expected = stats.chi2_contingency(ct) print(f"\nChi-square for {groups[0]} vs {groups[1]}") print(f"Chi2: {chi2:.2f}, p-value: {p:.4f}, DOF: {dof}") print("Contingency Table:\n", ct)

Run pairwise tests

pairwise_groups = ('Married', 'Never Married'), ('Married', 'Divorced'), ('Married', 'Cohabiting'), ('Never Married', 'Separated')

for g1, g2 in pairwise_groups: chi_square_pairwise(sub2, 'MARITALSTATUS', 'S2AQ8A', [g1, g2])

````````````````````````

Social anxiety is often seen as something people “grow out of” — maybe with age, confidence, or… education? It’s easy to assume that higher education, with all its exposure to social settings and professional environments, would reduce social anxiety. But what if that’s not actually the case?

I explored this exact question using the NESARC (National Epidemiologic Survey on Alcohol and Related Conditions) dataset. Although NESARC is primarily focused on alcohol and related conditions, it also provides validated and comprehensive measures of mental health disorders, including social anxiety.

🔍 The Test: Using ANOVA (Analysis of Variance) I compare social anxiety episodes across different education levels — from those who never finished high school to individuals with advanced degrees.

Here's the code: ----------------------------------------------------------------------

import pandas as pd import numpy as np import statsmodels.api as sm from statsmodels.formula.api import ols import matplotlib.pyplot as plt import seaborn as sns from statsmodels.stats.multicomp import pairwise_tukeyhsd

Load the data

data = pd.read_csv('nesarc.csv', low_memory=False) sub1 = data.copy()

Convert columns to numeric and coerce errors

sub1['S1Q6A'] = pd.to_numeric(sub1['S1Q6A'], errors='coerce') sub1['S7Q17C'] = pd.to_numeric(sub1['S7Q17C'], errors='coerce')

Drop rows with missing data in relevant columns

sub1 = sub1.dropna(subset=['S1Q6A', 'S7Q17C'])

Convert education to categorical

sub1['S1Q6A'] = sub1['S1Q6A'].astype('category')

Summary statistics

means = sub1.groupby('S1Q6A')['S7Q17C'].mean() stds = sub1.groupby('S1Q6A')['S7Q17C'].std() print("Group Means:\n", means) print("Group Standard Deviations:\n", stds)

ANOVA

model = ols('S7Q17C ~ C(S1Q6A)', data=sub1).fit() anova_table = sm.stats.anova_lm(model, typ=2) print("\nANOVA Table:\n", anova_table)

Tukey HSD for post hoc analysis

tukey = pairwise_tukeyhsd(endog=sub1['S7Q17C'], groups=sub1['S1Q6A'], alpha=0.05) print("\nTukey HSD Results:\n", tukey)

Boxplot for visualization

plt.figure(figsize=(12,6)) sns.boxplot(x='S1Q6A', y='S7Q17C', data=sub1) plt.xticks(rotation=90) plt.title('Social Anxiety Episodes by Education Level') plt.xlabel('Education Level') plt.ylabel('Number of Episodes') plt.tight_layout() plt.show() ---------------------------------------------------------------------------

📊 The Result: The p-value I found was 0.209 — meaning the differences observed between education levels were not statistically significant.

In plain terms: Your level of education doesn’t seem to predict how often you experience social anxiety. Whether you finished high school or earned a PhD, you're just as likely — or unlikely — to struggle with it.

And that’s an important insight.

🔎 Next Steps: Future studies could look into how employment status, social support, or therapy history intersect with social anxiety. The numbers might surprise us again.