Refraction (diopters) analysis
Juan Steibel
6/14/2021
Analysis of Refraction (diopters) datasets
library(tidyverse)
library(readxl)
library(lsmeans)
library(lme4)
library(ggplot2)
library(lmerTest)
setwd("C:/Users/marti/OneDrive/Documents/andras/AHFR/")Load data for Refraction, use graphics to summarize
Read data and show graphics
dts<-read_xlsx("Refraction and closest AScan for stats 5-4-2020.xlsx",sheet = "Sheet1")
table (dts$`Age at Refrection (days)`,dts$Status)##
## mutant wt
## 376 16 0
## 466 8 0
## 509 4 0
## 600 2 0
## 656 2 0
## 688 0 2
## 700 0 4
## 908 0 2
## 931 0 1
## 985 0 1
## 1067 0 4
## 1309 3 0
## 1535 0 1
## 1541 0 1
## 1591 2 0
## 1661 1 0
## 2188 2 0table (dts$`Age at Refraction (months)`,dts$Status)##
## mutant wt
## 12.3616438356164 16 0
## 15.3205479452055 8 0
## 16.7342465753425 4 0
## 19.7260273972603 2 0
## 21.5671232876712 2 0
## 22.6191780821918 0 2
## 23.013698630137 0 4
## 29.8520547945205 0 2
## 30.6082191780822 0 1
## 32.3835616438356 0 1
## 35.0794520547945 0 4
## 43.0356164383562 3 0
## 50.4657534246575 0 1
## 50.6630136986301 0 1
## 52.3068493150685 2 0
## 54.6082191780822 1 0
## 71.9342465753425 2 0ggplot(dts,aes(y=`Refraction (diopters)`,x=`Age at Ascan (days)`,color=`Animal ID`))+geom_point()+facet_grid(rows=vars(`Eye`),cols=vars(`Status`))
ggplot(dts,aes(y=`Refraction (diopters)`,x=`Age at Ascan (days)`))+geom_point()+facet_grid(rows=vars(`Eye`),cols=vars(`Status`))+geom_smooth(method="lm")## `geom_smooth()` using formula 'y ~ x'
I some evidence of age trend, to play safe let’s include age as a covariate and also its interaction with group.
ANOVA and mean comparisons
Similar to our previous analyses, a linear mixed model that has included fixed effects of status, eye, age and age*group term plus a random term for dog to account for repeated measures.
dts<-mutate(dts,age=as.numeric(`Age at Ascan (days)`)-mean(`Age at Ascan (days)`))
ah_lm<-lmer(`Refraction (diopters)`~Status+`Eye`+age*Status+(1|`Animal ID`),data=dts)
print(ah_lm)## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: `Refraction (diopters)` ~ Status + Eye + age * Status + (1 |
## `Animal ID`)
## Data: dts
## REML criterion at convergence: 215.8616
## Random effects:
## Groups Name Std.Dev.
## Animal ID (Intercept) 0.6292
## Residual 1.2344
## Number of obs: 56, groups: Animal ID, 32
## Fixed Effects:
## (Intercept) Statuswt EyeOS age Statuswt:age
## -5.327546 5.509551 -0.041093 0.000745 -0.002459anova(ah_lm,ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Status 175.406 175.406 1 26.184 115.1064 4.42e-11 ***
## Eye 0.022 0.022 1 27.161 0.0145 0.90511
## age 0.735 0.735 1 38.193 0.4826 0.49146
## Status:age 4.683 4.683 1 38.442 3.0732 0.08756 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lsmeans(ah_lm,pairwise~Status)## NOTE: Results may be misleading due to involvement in interactions## $lsmeans
## Status lsmean SE df lower.CL upper.CL
## mutant -5.348 0.241 26.8 -5.842 -4.85
## wt 0.161 0.454 26.0 -0.771 1.09
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mutant - wt -5.51 0.514 26.2 -10.729 <.0001
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-rogerCorrelation analyses
How are A-scan measures related to refraction measures?
cr<-cor(dts[,10:14])
cr[1,]## Refraction (diopters) AXL (mm) ACD (mm)
## 1.00000000 -0.24973978 -0.11523191
## LT (mm) VCD (mm)
## -0.20930949 -0.08740691small correlations.
dt<-dts[,10:14]
lmm<-lm(dt$`Refraction (diopters)`~as.matrix(dt[,2:5]))
summary(lmm)##
## Call:
## lm(formula = dt$`Refraction (diopters)` ~ as.matrix(dt[, 2:5]))
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.1516 -1.6335 -0.7898 2.5117 5.2030
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.75 16.76 1.655 0.104
## as.matrix(dt[, 2:5])AXL (mm) 17.72 52.37 0.338 0.736
## as.matrix(dt[, 2:5])ACD (mm) -18.67 52.58 -0.355 0.724
## as.matrix(dt[, 2:5])LT (mm) -19.92 52.65 -0.378 0.707
## as.matrix(dt[, 2:5])VCD (mm) -19.02 52.54 -0.362 0.719
##
## Residual standard error: 2.773 on 51 degrees of freedom
## Multiple R-squared: 0.0746, Adjusted R-squared: 0.002019
## F-statistic: 1.028 on 4 and 51 DF, p-value: 0.4019pairs(dt)
This does not help either. One last thing that could be atempted is some sort of non-linear relation. But without any prior knowledge, it would be too easy to overfit the data.
We can discuss this on THU.
Analysis of A-scan measures
print("AXL (mm)")## [1] "AXL (mm)"ah_lm<-lmer(`AXL (mm)`~Status+`Eye`+age*Status+(1|`Animal ID`),data=dts)
print(ah_lm)## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: `AXL (mm)` ~ Status + Eye + age * Status + (1 | `Animal ID`)
## Data: dts
## REML criterion at convergence: 102.4272
## Random effects:
## Groups Name Std.Dev.
## Animal ID (Intercept) 0.4021
## Residual 0.3101
## Number of obs: 56, groups: Animal ID, 32
## Fixed Effects:
## (Intercept) Statuswt EyeOS age Statuswt:age
## 20.7867158 -0.4737377 0.1838459 0.0003184 -0.0003355anova(ah_lm,ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Status 0.47400 0.47400 1 27.312 4.9293 0.03489 *
## Eye 0.42037 0.42037 1 25.054 4.3715 0.04685 *
## age 0.02988 0.02988 1 32.441 0.3107 0.58109
## Status:age 0.03681 0.03681 1 32.706 0.3828 0.54038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lsmeans(ah_lm,pairwise~Status)## NOTE: Results may be misleading due to involvement in interactions## $lsmeans
## Status lsmean SE df lower.CL upper.CL
## mutant 20.9 0.0997 27.5 20.7 21.1
## wt 20.4 0.1887 27.3 20.0 20.8
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mutant - wt 0.474 0.213 27.3 2.220 0.0349
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-rogerprint("ACD (mm)")## [1] "ACD (mm)"ah_lm<-lmer(`ACD (mm)`~Status+`Eye`+age*Status+(1|`Animal ID`),data=dts)
print(ah_lm)## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: `ACD (mm)` ~ Status + Eye + age * Status + (1 | `Animal ID`)
## Data: dts
## REML criterion at convergence: 97.5507
## Random effects:
## Groups Name Std.Dev.
## Animal ID (Intercept) 0.3633
## Residual 0.3061
## Number of obs: 56, groups: Animal ID, 32
## Fixed Effects:
## (Intercept) Statuswt EyeOS age Statuswt:age
## 4.348e+00 -2.155e-01 2.048e-01 -1.766e-05 -4.146e-05anova(ah_lm,ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Status 0.11164 0.11164 1 27.217 1.1913 0.28463
## Eye 0.52389 0.52389 1 25.277 5.5904 0.02604 *
## age 0.00218 0.00218 1 32.997 0.0233 0.87970
## Status:age 0.00063 0.00063 1 33.277 0.0067 0.93507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lsmeans(ah_lm,pairwise~Status)## NOTE: Results may be misleading due to involvement in interactions## $lsmeans
## Status lsmean SE df lower.CL upper.CL
## mutant 4.45 0.0923 27.5 4.26 4.64
## wt 4.24 0.1747 27.2 3.88 4.59
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mutant - wt 0.216 0.197 27.2 1.091 0.2846
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-rogerprint("LT (mm)")## [1] "LT (mm)"ah_lm<-lmer(`ACD (mm)`~Status+`Eye`+age*Status+(1|`Animal ID`),data=dts)
print(ah_lm)## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: `ACD (mm)` ~ Status + Eye + age * Status + (1 | `Animal ID`)
## Data: dts
## REML criterion at convergence: 97.5507
## Random effects:
## Groups Name Std.Dev.
## Animal ID (Intercept) 0.3633
## Residual 0.3061
## Number of obs: 56, groups: Animal ID, 32
## Fixed Effects:
## (Intercept) Statuswt EyeOS age Statuswt:age
## 4.348e+00 -2.155e-01 2.048e-01 -1.766e-05 -4.146e-05anova(ah_lm,ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Status 0.11164 0.11164 1 27.217 1.1913 0.28463
## Eye 0.52389 0.52389 1 25.277 5.5904 0.02604 *
## age 0.00218 0.00218 1 32.997 0.0233 0.87970
## Status:age 0.00063 0.00063 1 33.277 0.0067 0.93507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lsmeans(ah_lm,pairwise~Status)## NOTE: Results may be misleading due to involvement in interactions## $lsmeans
## Status lsmean SE df lower.CL upper.CL
## mutant 4.45 0.0923 27.5 4.26 4.64
## wt 4.24 0.1747 27.2 3.88 4.59
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mutant - wt 0.216 0.197 27.2 1.091 0.2846
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-rogerprint("VCD (mm)")## [1] "VCD (mm)"ah_lm<-lmer(`ACD (mm)`~Status+`Eye`+age*Status+(1|`Animal ID`),data=dts)
print(ah_lm)## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: `ACD (mm)` ~ Status + Eye + age * Status + (1 | `Animal ID`)
## Data: dts
## REML criterion at convergence: 97.5507
## Random effects:
## Groups Name Std.Dev.
## Animal ID (Intercept) 0.3633
## Residual 0.3061
## Number of obs: 56, groups: Animal ID, 32
## Fixed Effects:
## (Intercept) Statuswt EyeOS age Statuswt:age
## 4.348e+00 -2.155e-01 2.048e-01 -1.766e-05 -4.146e-05anova(ah_lm,ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Status 0.11164 0.11164 1 27.217 1.1913 0.28463
## Eye 0.52389 0.52389 1 25.277 5.5904 0.02604 *
## age 0.00218 0.00218 1 32.997 0.0233 0.87970
## Status:age 0.00063 0.00063 1 33.277 0.0067 0.93507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lsmeans(ah_lm,pairwise~Status)## NOTE: Results may be misleading due to involvement in interactions## $lsmeans
## Status lsmean SE df lower.CL upper.CL
## mutant 4.45 0.0923 27.5 4.26 4.64
## wt 4.24 0.1747 27.2 3.88 4.59
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mutant - wt 0.216 0.197 27.2 1.091 0.2846
##
## Results are averaged over the levels of: Eye
## Degrees-of-freedom method: kenward-roger