Analysis of Corneal Thickness
Juan Steibel
4/29/2020
Setup packages, please skip to the next section
## -- Attaching packages ------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.0 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.5
## v tidyr 1.0.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ---------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: emmeans## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepDiagnostics and graphics
Graphical output per dog indicates that the linear trend vary by dog. some have increasing trends, some decreasing trends.
There are 17 to 30 measures per dog and it’s always balanced between treatments
DPI varies between -8 and 1070. We should decide what to do with the negative dpi.
DPI per dog reveals that many dogs don’t have a 0 DPI and they don’t reach to 720, so the timepoints should be reviewed.
setwd("C:/Users/marti/OneDrive/Documents/andras/Kristin")
dts<-read_xlsx("All Dogs_Pachymetry_one sheet_STATS 2020.xlsx",sheet = "Pachy_all_dogs")
dts## # A tibble: 318 x 5
## Dog `psuedo title` eye `days post injection` µm
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Epoxy Dog 1 treated 29 583
## 2 Epoxy Dog 1 treated 43 582
## 3 Epoxy Dog 1 treated 63 575
## 4 Epoxy Dog 1 treated 126 548
## 5 Epoxy Dog 1 treated 154 532
## 6 Epoxy Dog 1 treated 167 543
## 7 Epoxy Dog 1 treated 196 532
## 8 Epoxy Dog 1 treated 224 591
## 9 Epoxy Dog 1 treated 245 561
## 10 Epoxy Dog 1 treated 266 547
## # ... with 308 more rowscolnames(dts)<-c("ID","dog","trt","dpi","CT")
dts## # A tibble: 318 x 5
## ID dog trt dpi CT
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Epoxy Dog 1 treated 29 583
## 2 Epoxy Dog 1 treated 43 582
## 3 Epoxy Dog 1 treated 63 575
## 4 Epoxy Dog 1 treated 126 548
## 5 Epoxy Dog 1 treated 154 532
## 6 Epoxy Dog 1 treated 167 543
## 7 Epoxy Dog 1 treated 196 532
## 8 Epoxy Dog 1 treated 224 591
## 9 Epoxy Dog 1 treated 245 561
## 10 Epoxy Dog 1 treated 266 547
## # ... with 308 more rowstable(dts$ID,dts$trt)##
## treated untreated
## BC 17 17
## Element 24 24
## Emulsion 30 30
## Epoxy 27 27
## Gaston 20 20
## Goldi 21 21
## Granny 20 20summary(dts$dpi)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.0 160.0 273.0 339.2 524.8 1071.0by(dts$dpi,dts$ID,summary)## dts$ID: BC
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.0 160.0 243.0 293.4 498.0 642.0
## ------------------------------------------------------------
## dts$ID: Element
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.0 140.8 244.0 339.8 498.0 1042.0
## ------------------------------------------------------------
## dts$ID: Emulsion
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.0 167.0 347.0 411.1 610.0 1071.0
## ------------------------------------------------------------
## dts$ID: Epoxy
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 29.0 203.0 323.0 366.8 557.2 730.0
## ------------------------------------------------------------
## dts$ID: Gaston
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.0 155.0 257.5 306.8 504.5 671.0
## ------------------------------------------------------------
## dts$ID: Goldi
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.0 146.0 243.0 299.1 498.0 671.0
## ------------------------------------------------------------
## dts$ID: Granny
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.0 155.0 258.0 307.0 504.5 671.0ggplot(dts,aes(y=CT,x=dpi,color=trt))+facet_wrap(.~ID)+geom_point()+geom_smooth(method = "lm")## `geom_smooth()` using formula 'y ~ x'
ggplot(dts,aes(y=CT,x=dpi,color=trt))+geom_point()+geom_smooth(method = "lm")## `geom_smooth()` using formula 'y ~ x'
Model selection
First we compare three models and the model with dog-specific intercept and slope is the best fit.
This model accounts for within-dog heteroneous variances and correlations over time.
dts<-mutate(dts,time=as.factor(dpi),x=(dpi-min(dpi))/(max(dpi)-min(dpi)))
m1<-lmer(CT~trt+dpi+trt*dpi+(1|ID),data=dts)
m2<-lmer(CT~trt+dpi+trt*dpi+(1|ID)+(x-1|ID),data=dts)
m3<-lmer(CT~trt+dpi+trt*dpi+(1+x|ID),data=dts)
anova(m1,m2)## refitting model(s) with ML (instead of REML)## Data: dts
## Models:
## m1: CT ~ trt + dpi + trt * dpi + (1 | ID)
## m2: CT ~ trt + dpi + trt * dpi + (1 | ID) + (x - 1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 6 2929.9 2952.4 -1458.9 2917.9
## m2 7 2887.7 2914.0 -1436.8 2873.7 44.164 1 3.019e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m3,m2)## refitting model(s) with ML (instead of REML)## Data: dts
## Models:
## m2: CT ~ trt + dpi + trt * dpi + (1 | ID) + (x - 1 | ID)
## m3: CT ~ trt + dpi + trt * dpi + (1 + x | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 7 2887.7 2914.0 -1436.8 2873.7
## m3 8 2888.0 2918.1 -1436.0 2872.0 1.7008 1 0.1922Analysis of significance
No interaction or main effect is significant in this analysis. We should not look into time-specific effects.
resu2<-m3
anova(resu2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## trt 1106.30 1106.30 1 302.534 2.6045 0.1076
## dpi 773.93 773.93 1 6.165 1.8220 0.2245
## trt:dpi 251.16 251.16 1 302.534 0.5913 0.4425lsmeans(resu2,"trt",at=list(dpi=0))## NOTE: Results may be misleading due to involvement in interactions## trt lsmean SE df lower.CL upper.CL
## treated 577 14.1 6.25 543 611
## untreated 583 14.1 6.25 549 618
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95lsmeans(resu2,"trt",at=list(dpi=120))## NOTE: Results may be misleading due to involvement in interactions## trt lsmean SE df lower.CL upper.CL
## treated 580 15 6.13 543 617
## untreated 585 15 6.13 549 622
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95lsmeans(resu2,"trt",at=list(dpi=360))## NOTE: Results may be misleading due to involvement in interactions## trt lsmean SE df lower.CL upper.CL
## treated 586 17.4 6.05 544 629
## untreated 590 17.4 6.05 547 632
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95lsmeans(resu2,"trt",at=list(dpi=720))## NOTE: Results may be misleading due to involvement in interactions## trt lsmean SE df lower.CL upper.CL
## treated 595 21.9 6.12 542 649
## untreated 596 21.9 6.12 543 650
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95