import pandas as pd

oasis_df = pd.read_csv('oasis_longitudinal.csv')

oasis_df = oasis_df.drop(['MRI ID','MR Delay', 'Hand'], axis=1)
oasis_df = oasis_df.rename({'Subject ID': 'Subject_ID', "M/F" : 'Gender'}, axis=1)
oasis_df['Gender'] = oasis_df['Gender'].astype('category')
oasis_df['Gender'] = oasis_df['Gender'].cat.codes

oasis_df

one_visit = oasis_df[oasis_df['Visit'] == 1] # subset of data
dummies = pd.get_dummies(one_visit["Group"]) # dummy coding for group
one_visit = one_visit.merge(dummies, left_index=True, right_index=True).dropna()
one_visit = one_visit[ one_visit["Demented"] != 1 ] # removed participants who already have dementia
one_visit.head()

both_visits = oasis_df[oasis_df['Visit'] <= 2].dropna()

both_visits = both_visits.groupby('Subject_ID').filter(lambda x: len(x) == 2)
both_visits_EDUC_SES_gender = both_visits[['Subject_ID','EDUC','SES','Gender']] #keep values that do not change
both_visits = both_visits.sort_values(["Subject_ID", "Visit"])
both_visits_diff = both_visits.groupby(['Subject_ID'])[['Age', 'MMSE',
                                                        'CDR', 'eTIV', 'nWBV', 'ASF']].diff().dropna()
both_visits = both_visits[["Subject_ID", "Group"]].merge(both_visits_diff, left_index=True, right_index=True)
both_visits = both_visits.merge(both_visits_EDUC_SES_gender, left_index=True, right_index=True, suffixes=(None,'_y'))
both_visits = both_visits.drop(columns = ['Subject_ID_y'])

dummies = pd.get_dummies(both_visits["Group"])
both_visits = both_visits.merge(dummies, left_index=True, right_index=True).dropna()

both_visits = both_visits[both_visits["Demented"] != 1]
both_visits = both_visits.drop(columns='Group')

import numpy as np
import statsmodels.formula.api as smf
import statsmodels.api as sm
import pylab as py
import seaborn as sns
import matplotlib.pyplot as plt

both_visits.head()

model_fit = smf.ols('Converted ~ Age+EDUC+SES+MMSE+eTIV+nWBV+ASF+Gender', data=both_visits).fit()

X = pd.DataFrame(both_visits, columns=['Age','EDUC','SES', 'MMSE','eTIV','nWBV','ASF', 'Gender'])
y = pd.DataFrame(both_visits.Converted)
dataframe = pd.concat([X, y], axis=1)

table = sm.stats.anova_lm(model_fit, typ=2)
with pd.option_context('display.max_rows', None, 'display.max_columns', None):  
    display(table.sort_values("F", ascending=False))

#model values
model_fitted_y = model_fit.fittedvalues
model_residuals = model_fit.resid
model_norm_residuals = model_fit.get_influence().resid_studentized_internal
model_norm_residuals_abs_sqrt = np.sqrt(np.abs(model_norm_residuals))
model_abs_resid = np.abs(model_residuals)
model_leverage = model_fit.get_influence().hat_matrix_diag
model_cooks = model_fit.get_influence().cooks_distance[0]

from statsmodels.graphics.gofplots import ProbPlot

QQ = ProbPlot(model_norm_residuals)
plot_lm_2 = QQ.qqplot(line='45', alpha=0.5, color='#4C72B0', lw=1)
plot_lm_2.axes[0].set_title('Normal Q-Q Plot')
plot_lm_2.axes[0].set_xlabel('Theoretical Quantiles')
plot_lm_2.axes[0].set_ylabel('Standardized Residuals');

from scipy.stats import shapiro
from tabulate import tabulate

shapiro_results = []

for column in both_visits:
    try:
        s = shapiro(both_visits[column])
        shapiro_results.append([column, s.statistic, s.pvalue])
    except Exception as e:
        pass

df = pd.DataFrame(shapiro_results)
df = df.drop([9,10,11], axis=0) #drop dummy measures
print('Shapiro-Wilk Test \n')
print(tabulate(df, headers=['Measure', 'W statistic', 'p value'],showindex="never",tablefmt="presto"))  

#fit model
mixed = smf.mixedlm('Converted ~ MMSE*nWBV*eTIV*EDUC*SES',data=both_visits, groups=both_visits['Subject_ID'])
mixed_fit = mixed.fit(reml=False)

#print summary
mixed_fit_results = mixed_fit.summary()
print(mixed_fit_results)
BIC_1 = mixed_fit.bic
print('Bayesian Information Criteria (BIC): ', BIC_1)

#convert summary results to dataframe
def summary_to_df(results):
    pvals = results.pvalues
    coeff = results.params
    conf_lower = results.conf_int()[0] #lower confidence interval of the fitted parameters
    conf_higher = results.conf_int()[1]

    results_df = pd.DataFrame({"pvals":pvals,
                               "coeff":coeff,
                               "conf_lower":conf_lower,
                               "conf_higher":conf_higher
                                })

    results_df = results_df[["coeff","pvals","conf_lower","conf_higher"]]
    return results_df

summary_df = summary_to_df(mixed_fit)
sig_pvals = summary_df[summary_df ['pvals']<= 0.05]
sig_pvals

#fit model 2
mixed2 = smf.mixedlm('Converted ~ MMSE+SES+EDUC+nWBV+eTIV+Age+Gender',data=both_visits, groups=both_visits['Subject_ID'])
mixed2_fit = mixed2.fit(reml=False)

print(mixed2_fit.summary())
print('Bayesian Information Criteria (BIC): ', mixed2_fit.bic)

summary_df2 = summary_to_df(mixed2_fit)
sig_pvals2 = summary_df2[summary_df2['pvals'] <= 0.05]
sig_pvals2

oasis_df2 = oasis_df.dropna() #remove missing scores

g = sns.regplot(x='SES',y='EDUC', data=one_visit,color="skyblue", scatter_kws={'s':30}, label = 'A');
g = g.set(xlim=(0,6), xlabel="SES (1-5)", ylabel="Education (Years)")


correlation_matrix = np.corrcoef(one_visit['SES'], one_visit['EDUC'])
correlation_xy = correlation_matrix[0,1]
r_squared = correlation_xy**2


print("R-Squared:","{:.2f}".format(r_squared))

g2 = sns.regplot(x='SES',y='MMSE', data=one_visit);
g2 = g2.set(xlim=(0, 6))

correlation_matrix = np.corrcoef(one_visit['SES'], one_visit['MMSE'])
correlation_xy = correlation_matrix[0,1]
r_squared = correlation_xy**2

print("R-Squared:","{:.4f}".format(r_squared))

g3 = sns.regplot(x='EDUC',y='MMSE', data=one_visit, color="g");
g3 = g3.set(xlim=(7, 25))

oasis_df2 = one_visit.dropna()
correlation_matrix = np.corrcoef(one_visit['EDUC'], one_visit['MMSE'])
correlation_xy = correlation_matrix[0,1]
r_squared = correlation_xy**2

print("R-Squared:","{:.2f}".format(r_squared))