Examining the data for 5 ADHD participants (AD01, AD04, AD08, AD10, and AD15) before and after medication

This section examines whether taking ADHD medicine before or after —- has an effect on the response variables while controlling for other factors.

# data can be easily uploaded via the 'import dataset' tab in R studio...
library(readxl)
adhd.med <- read_excel(file.choose(), sheet = "ADHDmed")


adhd.med2 <- subset(adhd.med, select=c(
# independent variables
  `ID entered in RETeval device`,age,gender,VERT,
  Ethnicity_Simple,`Tested eye`,iris,
  Strength,Before_After,
# dependent variables
  a_time,   a_amp, b_time,  b_amp, bdiva_amp,   p72, Tmin, BT,
  p_ratio,  w_ratio))
  
# checking the data's structure...
str(adhd.med2)
## tibble [355 x 19] (S3: tbl_df/tbl/data.frame)
##  $ ID entered in RETeval device: chr [1:355] "AD01" "AD01" "AD01" "AD01" ...
##  $ age                         : num [1:355] 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 ...
##  $ gender                      : num [1:355] 0 0 0 0 0 0 0 0 0 0 ...
##  $ VERT                        : num [1:355] 2 2 2 2 2 2 2 2 2 2 ...
##  $ Ethnicity_Simple            : num [1:355] 4 4 4 4 4 4 4 4 4 4 ...
##  $ Tested eye                  : num [1:355] 1 1 1 1 1 1 1 1 1 1 ...
##  $ iris                        : num [1:355] 1.38 1.38 1.38 1.38 1.38 ...
##  $ Strength                    : num [1:355] -0.367 -0.367 -0.119 0.114 0.398 0.398 0.477 0.477 0.602 0.602 ...
##  $ Before_After                : num [1:355] 0 0 0 0 0 0 0 0 0 0 ...
##  $ a_time                      : num [1:355] 14.9 15.1 13.3 11.3 12.6 ...
##  $ a_amp                       : num [1:355] -4.22 -3.01 -3.3 -5.79 -7.2 ...
##  $ b_time                      : num [1:355] 20.2 21.8 21.7 24.2 24.5 ...
##  $ b_amp                       : num [1:355] 14.2 10.7 17.2 25.9 38.1 ...
##  $ bdiva_amp                   : num [1:355] 3.37 3.56 5.2 4.48 5.29 ...
##  $ p72                         : num [1:355] 1.74 -1.02 -2.67 -7.69 -4.42 ...
##  $ Tmin                        : num [1:355] 99.9 99.4 100.1 100.1 85.9 ...
##  $ BT                          : num [1:355] -2.63 -6.23 -4.41 -11.78 -8.71 ...
##  $ p_ratio                     : num [1:355] -0.175 0.133 0.192 0.382 0.143 ...
##  $ w_ratio                     : num [1:355] 0.888 1.301 1.064 1.231 1.04 ...
dim(adhd.med2)
## [1] 355  19
# rendering variables into their right type
adhd.med2$`ID entered in RETeval device` <- as.factor(adhd.med2$`ID entered in RETeval device`)
adhd.med2$gender <- as.factor(adhd.med2$gender)
adhd.med2$VERT <- as.factor(adhd.med2$VERT)
   adhd.med2$VERT<-relevel(adhd.med2$VERT, ref='2') # ref level in this variable
adhd.med2$Ethnicity_Simple <- as.factor(adhd.med2$Ethnicity_Simple)
adhd.med2$`Tested eye` <- as.factor(adhd.med2$`Tested eye`)
adhd.med2$Strength <- as.factor(adhd.med2$Strength)
adhd.med2$Before_After <- as.factor(adhd.med2$Before_After)


# changing the name of the columns
colnames(adhd.med2) <- c('Participant','Age','Gender','Vert',
                'Ethnicity','Eye','Iris','Strength','med before-after',
                'a_time', 'a_amp',  'b_time',   'b_amp','bdiva_amp',
                'p72',  'Tmin', 'BT',   'p_ratio',  'w_ratio')


# relabeling some of the variables' levels
library(plyr)
adhd.med2$Gender <- 
mapvalues(adhd.med2$Gender, 
  from = c('0','1'), 
  to = c('Male','Female'))

adhd.med2$Ethnicity <- 
mapvalues(adhd.med2$Ethnicity, 
  from = c('1','2','3','4','5'), 
  to = c('Caucasian','Asian','Afro-Caribbean','Latino','Mixed'))

adhd.med2$Eye <-
mapvalues(adhd.med2$Eye, 
  from = c('0','1'), 
  to = c('Right','Left'))

adhd.med2$`med before-after` <-
mapvalues(adhd.med2$`med before-after`, 
  from = c('0','1'), 
  to = c('Before','After'))


# checking if there are missing values...
# see https://cran.r-project.org/web/packages/naniar/vignettes/naniar-visualisation.html
# for more visualisations
#library(naniar)
#gg_miss_upset(adhd.med2)

Yep definitely a good thing to not have NAs!.

Descriptive statistics

This section shows some descriptive statistics along with robust estimations of location (median) and scale (MAD = median absolute deviations) for all dependent variables regardless of Strength and Eye.

library(tidyverse)

# Note 'Eye' is left out in order to avoid duplicates
# VERT and Strength are also not included
some_stats_by_participants.adhd <- adhd.med2 %>% 
  group_by(Participant, Age, Gender, Ethnicity, `med before-after`) %>% 
  summarize(median.A.time = median(a_time),
            MAD.A.time=mad(a_time),
            CV.A.time = mad(a_time)/median(a_time),
            
            median.A.amp = median(a_amp),
            MAD.A.amp=mad(a_amp),
            CV.A.amp = mad(a_amp)/median(a_amp),
            
            median.B.time = median(b_time),
            MAD.B.time=mad(b_time),
            CV.B.time = mad(b_time)/median(b_time),
            
            median.B.amp = median(b_amp),
            MAD.B.amp=mad(b_amp),
            CV.B.amp = mad(b_amp)/median(b_amp),
            
            median.bdiva_amp = median(bdiva_amp),
            MAD.bdiva_amp=mad(bdiva_amp),
            CV.bdiva_amp = mad(bdiva_amp)/median(bdiva_amp),
            
            median.p72 = median(p72),
            MAD.p72 =mad(p72),
            CV.p72 = mad(p72)/median(p72),
            
            median.Tmin = median(Tmin),
            MAD.Tmin =mad(Tmin),
            CV.Tmin = mad(Tmin)/median(Tmin),
            
             median.BT = median(BT),
            MAD.BT=mad(BT),
            CV.BT = mad(BT)/median(BT),
            
            median.p_ratio = median(p_ratio),
            MAD.p_ratio =mad(p_ratio),
            CV.p_ratio = mad(p_ratio)/median(p_ratio),
            
            median.w_ratio = median(w_ratio),
            MAD.w_ratio =mad(w_ratio),
            CV.w_ratio = mad(w_ratio)/median(w_ratio))

# filtering the data for demographics only
demographics.adhd <- some_stats_by_participants.adhd[ , c(1:5)]


# number of observations per gender, ethnicity and group
# note: observations can come from the same participant
library(ggplot2)
library(ggpubr)
theme_set(theme_pubr())

count.data.adhd<-with(demographics.adhd, 
                      ftable(Gender, `med before-after`, Ethnicity))

count.data.adhd1<-as.data.frame(count.data.adhd) 
names(count.data.adhd1)[4] <- c('Count')

As in the above analyses, an important aspect here is that there should be one iris measure per eye per participant. That is, each participant should have two iris measurements; one for the left eye and another for the right eye.

# here b-time is simply an excuse to obtain unique values for the other variables
# the way to do is by setting the function to 'class'
iris.aggregate.adhd<-aggregate(data=adhd.med2, 
    b_time ~ Participant + Gender + Age + Ethnicity + `med before-after` + Eye + Iris, FUN='class')

iris.aggregate.sorted.adhd<-
  data.frame(sort(table(iris.aggregate.adhd$Participant)))

# those with more than four observations
iris.aggregate.sorted.adhd$Var1[iris.aggregate.sorted.adhd$Freq>4]
## factor(0)
## Levels: AD01 AD04 AD08 AD10 AD15
# those with less than 2 observations
iris.aggregate.sorted.adhd$Var1[iris.aggregate.sorted.adhd$Freq<4]
## [1] AD01
## Levels: AD01 AD04 AD08 AD10 AD15

Modelling the effect of ‘before’ vs ‘after’ medicine

The effect of ‘before’ vs ‘after’ medicine in five ADHD participants will be assessed on each dependent variable while controlling for other variables.

a-time

A GEE model gives error. According to some sites it’s due to the design matrix not being invertible and this happens due to high multicolinearity or from the \(n < p\) issue; i.e. less observations than variables 1. Recall, the data is from just 5 participants.

A robust linear fixed effects models will be used 2.

library(robustbase)
library(knitr)

a.time.model.adhd <- lmrob(a_time ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

library(car)
anova.a.time.model.adhd <- Anova(a.time.model.adhd,
                                 test.statistic='Chisq', type=3)
anova.a.time.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 4.4418294 0.0350687
med before-after 1 0.0501923 0.8227290
Strength 9 136.0113643 0.0000000
Vert 2 12.5818053 0.0018531
Ethnicity 2 4.4761873 0.1066616
Iris 1 0.9164068 0.3384202
Gender 1 0.0385071 0.8444283
Age 1 0.1029108 0.7483641
Eye 1 3.5289429 0.0603062
med before-after:Strength 9 2.8845030 0.9686984

a-amp

a.amp.model.adhd <- lmrob(a_amp ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.a.amp.model.adhd <- Anova(a.amp.model.adhd,
                                test.statistic='Chisq', type=3)
anova.a.amp.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 13.2633800 0.0002706
med before-after 1 0.1431994 0.7051209
Strength 9 170.9715197 0.0000000
Vert 2 10.0895148 0.0064430
Ethnicity 2 12.8083152 0.0016547
Iris 1 0.8250441 0.3637095
Gender 1 23.0495922 0.0000016
Age 1 7.8366010 0.0051199
Eye 1 1.4759473 0.2244097
med before-after:Strength 9 5.8626030 0.7535870

b-time

b.time.model.adhd <- lmrob(b_time ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.b.time.model.adhd <- Anova(b.time.model.adhd,
                                 test.statistic='Chisq', type=3)
anova.b.time.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 46.4339618 0.0000000
med before-after 1 1.7290445 0.1885330
Strength 9 4785.2874645 0.0000000
Vert 2 1.6592069 0.4362222
Ethnicity 2 91.9502639 0.0000000
Iris 1 0.0108317 0.9171096
Gender 1 6.0001896 0.0143043
Age 1 5.3688885 0.0204990
Eye 1 2.4566369 0.1170291
med before-after:Strength 9 5.6545118 0.7739360

b-amp

b.amp.model.adhd <- lmrob(b_amp ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.b.amp.model.adhd <- Anova(b.amp.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.b.amp.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 89.0104283 0.0000000
med before-after 1 0.7861413 0.3752695
Strength 9 672.2421771 0.0000000
Vert 2 14.6578668 0.0006563
Ethnicity 2 113.2517539 0.0000000
Iris 1 0.0001358 0.9907008
Gender 1 155.6425186 0.0000000
Age 1 87.6112895 0.0000000
Eye 1 7.3469220 0.0067178
med before-after:Strength 9 12.9368234 0.1654858

bdiva_amp

bdiva_amp.model.adhd <- lmrob(bdiva_amp ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.bdiva_amp.model.adhd <- Anova(bdiva_amp.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.bdiva_amp.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 0.2576293 0.6117536
med before-after 1 1.9772521 0.1596802
Strength 9 93.7781270 0.0000000
Vert 2 2.5760579 0.2758139
Ethnicity 2 0.8432307 0.6559863
Iris 1 15.2468690 0.0000943
Gender 1 0.8547777 0.3552043
Age 1 2.0547066 0.1517365
Eye 1 0.5243572 0.4689889
med before-after:Strength 9 7.7402809 0.5605167

p72

p72.model.adhd <- lmrob(p72 ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.p72.model.adhd <- Anova(p72.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.p72.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 7.992669 0.0046967
med before-after 1 1.735135 0.1877564
Strength 9 77.338688 0.0000000
Vert 2 2.491695 0.2876970
Ethnicity 2 6.040445 0.0487904
Iris 1 5.524869 0.0187480
Gender 1 5.426438 0.0198341
Age 1 3.874992 0.0490105
Eye 1 1.655538 0.1982069
med before-after:Strength 9 9.859407 0.3619866

Tmin

Tmin.model.adhd <- lmrob(Tmin ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.Tmin.model.adhd <- Anova(Tmin.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.Tmin.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 9.2165863 0.0023983
med before-after 1 8.4136186 0.0037242
Strength 9 15.2563402 0.0841312
Vert 2 0.3297030 0.8480196
Ethnicity 2 116.1197838 0.0000000
Iris 1 0.7209381 0.3958364
Gender 1 8.2661757 0.0040390
Age 1 4.5640509 0.0326498
Eye 1 0.0003914 0.9842168
med before-after:Strength 9 9.6684205 0.3779891

BT

BT.model.adhd <- lmrob(BT ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.BT.model.adhd <- Anova(BT.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.BT.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 6.927651 0.0084873
med before-after 1 1.803102 0.1793380
Strength 9 52.678194 0.0000000
Vert 2 3.168554 0.2050960
Ethnicity 2 10.077884 0.0064806
Iris 1 6.754120 0.0093531
Gender 1 3.352360 0.0671087
Age 1 2.098274 0.1474655
Eye 1 3.746515 0.0529177
med before-after:Strength 9 12.310007 0.1963914

p_ratio

p_ratio.model.adhd <- lmrob(p_ratio ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.p_ratio.model.adhd <- Anova(p_ratio.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.p_ratio.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 0.2532536 0.6147934
med before-after 1 1.3677893 0.2421921
Strength 9 44.3435452 0.0000012
Vert 2 0.8132433 0.6658961
Ethnicity 2 0.8700534 0.6472474
Iris 1 1.6817016 0.1946985
Gender 1 0.6581762 0.4172045
Age 1 0.1000497 0.7517700
Eye 1 0.2108639 0.6460911
med before-after:Strength 9 10.8657120 0.2850262

w_ratio

w_ratio.model.adhd <- lmrob(w_ratio ~ 
              `med before-after` + 
    Strength + Vert + Ethnicity + Iris + Gender + Age + Eye
    + `med before-after`:Strength, 
                           data = adhd.med2)

anova.w_ratio.model.adhd <- Anova(w_ratio.model.adhd, 
                    test.statistic='Chisq', type=3)
anova.w_ratio.model.adhd %>% kable()
Df Chisq Pr(>Chisq)
(Intercept) 1 7.330515e+01 0.0000000
med before-after 1 3.492242e+03 0.0000000
Strength 9 4.715867e+04 0.0000000
Vert 2 8.735018e+00 0.0126828
Ethnicity 2 4.182045e+00 0.1235607
Iris 1 2.532768e+00 0.1115044
Gender 1 9.956382e-01 0.3183682
Age 1 7.701000e-04 0.9778614
Eye 1 5.023976e+00 0.0249987
med before-after:Strength 9 3.718770e+03 0.0000000

  1. https://stats.stackexchange.com/questions/76488/error-system-is-computationally-singular-when-running-a-glm??

  2. A good explanation in relation to the type of ANOVA used can be seen here: https://mcfromnz.wordpress.com/2011/03/02/anova-type-iiiiii-ss-explained/??