Andras Komaromy Lab - Vanessa Data analysis
JPS-V-A
July 18, 2019
Retrospective Data Analysis
Setup of datasets
library(tidyverse)## -- Attaching packages --------------------------------------------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.1.0 v purrr 0.3.2
## v tibble 2.1.1 v dplyr 0.8.0.1
## v tidyr 0.8.3 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## -- Conflicts ------------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(readxl)
library(ggplot2)
rm(list=ls())
ls()## character(0)setwd("C:/Users/marti/OneDrive/Documents/andras/Vanessa/")
dts<-read_xlsx("raphtis compiled IOP measurements.xlsx")Various checks
We fixed various errors in the labeling: Non-Affected= Normal? and ADAMTS10 carrier = ADAMTS10 Carrier. We also create new variables needed for analyses and continue to check data consistency.
dim(dts)## [1] 6955 12head(dts)## # A tibble: 6 x 12
## Dog Gender Status DOB Date `Age in Days`
## <chr> <chr> <chr> <dttm> <dttm> <dbl>
## 1 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## 2 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## 3 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## 4 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## 5 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## 6 Aldo M ADAMT~ 2012-03-02 00:00:00 2012-05-29 00:00:00 88
## # ... with 6 more variables: `Age in Months` <dbl>, Time <chr>,
## # Tonomoter <chr>, Eye <chr>, IOP <dbl>, Puppy <chr>table(dts$Dog)%>%length()## [1] 68table(dts$Dog)##
## Aldo Andy Angel Anja Anjie
## 132 96 90 144 18
## Anton Arlington Armin Arnold Baron
## 144 108 219 84 156
## Battle Creek Betty Britney Spears Capone Chubbie
## 54 150 6 204 102
## Cleopatra Czar Dandelion Daphne Dogwood
## 90 302 166 154 142
## Dryas Emulsion Esteban Fargo Felix
## 148 84 24 42 24
## Gaston Gepetto Gilda Goldilocks Granny
## 132 144 174 138 132
## Gretel Grimm Grumpy Harry Hermione
## 138 144 144 132 18
## Hope Hoss Howgwarts Javier Jax
## 108 36 138 24 24
## Jayne Kepler Kuiper Lacey Lady
## 102 58 34 36 186
## Laser Lathe Lepton Liberty Lumi
## 58 58 58 108 70
## Luminous Lunar Mabel Macho Mateo
## 18 18 28 198 24
## Mattawan Nadia Piper Porkchop Princess
## 60 114 180 36 66
## Red Poppy Rita Spinning Stefano Toby
## 108 120 126 24 162
## Turtle Dove Waterloo Yankee Doodle
## 180 108 108table(dts$Gender)##
## F M
## 3256 3699table(dts$Status)##
## ACHM Affected ADAMTS10 Affected ADAMTS10 Carrier Normal
## 1242 4111 618 144
## wt
## 840dts<-mutate(dts,Status=tolower(Status))
table(dts$Status) #fied##
## achm affected adamts10 affected adamts10 carrier normal
## 1242 4111 618 144
## wt
## 840#create a variable that labels Adamts10 status as: affected vs non_affected
dts<-mutate(dts,adams_affected=Status=="adamts10 affected")
table(dts$adams_affected)##
## FALSE TRUE
## 2844 4111table(dts$`Age in Days`)##
## 19 20 21 22 23 24 25 28 29 30 31 32 33 34 35
## 8 8 8 12 12 12 20 20 12 12 20 20 8 8 20
## 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## 20 12 12 20 8 50 20 12 12 20 20 8 8 20 12
## 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
## 12 20 20 8 8 20 12 18 20 20 14 8 20 12 12
## 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 20 20 8 8 20 12 12 20 20 14 8 20 6 12 20
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
## 20 8 50 38 12 12 20 56 8 8 20 12 12 20 24
## 96 97 98 99 103 109 110 111 112 113 115 116 117 118 123
## 8 8 20 18 6 36 14 24 6 18 42 24 6 6 24
## 124 125 126 128 130 131 132 137 141 142 150 154 155 162 169
## 6 6 18 6 36 6 6 48 6 42 6 24 12 6 18
## 171 172 173 175 180 181 182 183 196 197 201 206 207 208 209
## 42 24 12 6 6 24 36 6 8 6 6 42 24 12 24
## 217 224 225 231 237 238 248 251 252 253 258 259 264 266 269
## 11 6 6 24 42 6 6 6 6 6 12 36 6 6 42
## 271 272 274 279 280 281 282 283 284 285 287 288 289 290 291
## 6 24 6 12 24 6 6 6 6 12 30 6 6 6 6
## 292 293 294 295 296 297 298 299 300 301 302 304 305 308 309
## 18 6 24 12 12 54 6 6 12 36 18 6 6 24 6
## 312 313 314 315 318 319 322 324 325 326 327 329 334 335 336
## 6 6 6 24 42 12 18 6 6 18 6 6 6 12 6
## 337 340 341 342 343 350 356 357 358 360 361 363 364 365 372
## 48 6 6 6 6 12 6 6 6 12 24 6 6 6 6
## 377 378 379 383 384 386 390 393 394 396 397 402 403 406 410
## 18 36 48 6 24 4 6 6 6 6 12 12 6 6 30
## 411 415 417 425 431 432 435 437 438 443 451 453 456 458 459
## 12 36 12 6 12 12 6 18 24 6 48 6 6 6 18
## 460 462 463 465 466 469 470 473 474 475 477 478 479 484 487
## 6 12 6 6 12 12 6 6 18 24 10 12 42 6 12
## 488 491 495 497 499 500 501 503 505 510 512 513 514 516 519
## 6 6 6 6 6 6 6 6 6 6 6 12 54 6 6
## 522 523 524 526 528 529 533 536 537 540 542 544 549 550 556
## 6 6 18 6 6 6 6 6 6 6 42 6 6 6 6
## 558 560 563 566 567 568 569 570 574 577 581 586 588 593 594
## 6 12 6 12 12 6 6 6 6 42 6 6 12 12 18
## 595 598 603 604 605 609 610 611 615 616 621 622 623 626 630
## 12 6 42 6 6 6 12 6 12 12 6 12 6 42 12
## 634 636 637 641 645 646 647 649 652 655 657 658 663 664 670
## 6 6 6 6 12 6 6 6 42 6 12 10 6 12 6
## 677 679 682 683 684 685 686 693 698 705 707 709 710 711 715
## 6 1 6 6 6 6 6 42 6 6 6 12 6 6 12
## 719 720 721 723 730 733 736 740 741 744 745 747 751 761 766
## 12 48 6 6 12 6 6 6 18 6 30 6 6 6 24
## 768 769 773 774 776 777 779 796 800 803 804 808 809 811 818
## 6 12 6 6 6 6 6 6 6 6 18 6 6 6 6
## 822 825 829 830 834 835 839 845 852 858 867 870 871 878 879
## 6 6 6 6 6 6 12 6 6 6 6 6 6 6 6
## 883 884 899 900 908 919 924 925 927 939 946 948 951 952 960
## 6 6 6 6 6 18 12 6 6 6 6 6 6 6 6
## 963 967 971 972 973 976 985 992 994 996 997 998 1000 1011 1014
## 6 6 6 6 6 6 6 6 6 12 12 6 6 6 6
## 1017 1022 1039 1048 1049 1066 1070 1083 1085 1090 1091 1095 1096 1113 1116
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 1121 1125 1126 1132 1145 1148 1154 1157 1168 1174 1176 1180 1181 1185 1193
## 6 18 6 6 6 6 6 6 6 18 6 6 6 6 6
## 1195 1197 1202 1216 1220 1221 1223 1230 1243 1248 1258 1265 1269 1285 1288
## 24 6 6 6 6 6 18 12 6 6 24 12 6 6 6
## 1294 1295 1302 1311 1312 1317 1324 1330 1353 1358 1360 1366 1373 1380 1381
## 6 18 12 18 12 18 12 12 6 6 18 6 12 12 18
## 1393 1401 1402 1406 1408 1409 1421 1430 1431 1442 1444 1449 1455 1456 1465
## 6 18 6 6 12 18 6 6 6 18 18 12 6 6 6
## 1468 1475 1481 1490 1493 1497 1503 1525 1528 1531 1532 1539 1554 1557 1559
## 18 12 18 18 6 18 18 6 6 6 12 12 6 6 42
## 1560 1566 1572 1577 1587 1599 1600 1603 1606 1607 1616 1620 1627 1628 1644
## 6 24 6 6 24 6 12 6 6 12 6 6 6 18 6
## 1645 1653 1654 1655 1660 1671 1676 1678 1685 1690 1692 1702 1707 1709 1717
## 6 6 18 6 12 6 12 6 12 6 6 6 6 12 6
## 1718 1727 1744 1745 1749 1751 1772 1779 1786 1805 1807 1811 1814 1815 1819
## 18 6 12 30 6 12 6 12 18 6 6 6 12 6 3
## 1833 1834 1837 1839 1854 1857 1861 1863 1864 1874 1882 1888 1889 1891 1900
## 6 6 6 6 3 6 6 6 6 6 3 6 6 6 6
## 1912 1917 1922 1926 1930 1931 1943 1958 1962 1964 1966 1985 1991 1992 1993
## 6 3 6 6 6 6 3 6 12 6 3 6 12 3 6
## 1997 2008 2011 2020 2032 2033 2049 2060 2077 2085 2095 2099 2106 2112 2117
## 6 6 6 6 12 3 6 9 6 9 6 6 3 6 6
## 2132 2140 2154 2161 2176 2190 2197 2203 2210 2211 2225 2229 2238 2279 2281
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 2305 2323 2327 2330 2336 2344 2351 2369 2377 2396 2400 2437 2490 2491 2515
## 6 6 6 6 6 6 6 6 6 12 6 6 6 6 6
## 2539 2617
## 6 6#count records in dogs by adamsts10 status:
tb<-table(dts$Dog,dts$adams_affected)
print(tb)##
## FALSE TRUE
## Aldo 0 132
## Andy 0 96
## Angel 0 90
## Anja 0 144
## Anjie 18 0
## Anton 0 144
## Arlington 108 0
## Armin 0 219
## Arnold 0 84
## Baron 0 156
## Battle Creek 54 0
## Betty 0 150
## Britney Spears 6 0
## Capone 204 0
## Chubbie 0 102
## Cleopatra 90 0
## Czar 0 302
## Dandelion 0 166
## Daphne 0 154
## Dogwood 0 142
## Dryas 0 148
## Emulsion 0 84
## Esteban 24 0
## Fargo 0 42
## Felix 24 0
## Gaston 0 132
## Gepetto 0 144
## Gilda 174 0
## Goldilocks 0 138
## Granny 0 132
## Gretel 0 138
## Grimm 0 144
## Grumpy 0 144
## Harry 132 0
## Hermione 18 0
## Hope 108 0
## Hoss 0 36
## Howgwarts 138 0
## Javier 24 0
## Jax 24 0
## Jayne 102 0
## Kepler 0 58
## Kuiper 0 34
## Lacey 36 0
## Lady 186 0
## Laser 0 58
## Lathe 0 58
## Lepton 0 58
## Liberty 108 0
## Lumi 0 70
## Luminous 0 18
## Lunar 0 18
## Mabel 0 28
## Macho 198 0
## Mateo 24 0
## Mattawan 60 0
## Nadia 114 0
## Piper 180 0
## Porkchop 36 0
## Princess 0 66
## Red Poppy 108 0
## Rita 0 120
## Spinning 126 0
## Stefano 24 0
## Toby 0 162
## Turtle Dove 180 0
## Waterloo 108 0
## Yankee Doodle 108 0write.csv(tb,file="dog_by_status.csv")
status_dog=(tb>0)*1
status_dog[rowSums(status_dog)>1,]##
## FALSE TRUEstatus_dog[rowSums(status_dog)<1,]##
## FALSE TRUEcolSums(status_dog)## FALSE TRUE
## 31 37#check if age in months is computed correctly > OK (+/- 1 day!)
mutate(dts,months=`Age in Days`/30)%>%select(`Age in Months`,months)%>%mutate(diff=abs(`Age in Months`-months))%>%filter(diff>0.05)## # A tibble: 0 x 3
## # ... with 3 variables: `Age in Months` <dbl>, months <dbl>, diff <dbl>#far more AM observations due to puppies and some dogs not having been observed later in the day.
table(dts$Time)##
## AM NOON PM
## 3039 1965 1951table(dts$Tonomoter)##
## TonoPen TonoVet
## 1356 5599table(dts$Eye)##
## OD OS
## 3461 3494Descriptive analyses
summary(dts$IOP)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 4.00 13.00 16.00 16.60 19.33 67.00 2by(dts$IOP,dts$adams_affected,summary)## dts$adams_affected: FALSE
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 4.00 11.00 13.33 13.56 15.67 34.00 2
## --------------------------------------------------------
## dts$adams_affected: TRUE
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.0 15.0 18.0 18.7 22.0 67.0d1<-ggplot(data = dts,aes(x=IOP))+geom_density()
d2<-ggplot(data = dts,aes(x=IOP))+geom_histogram()+geom_vline(xintercept = 16.6)
d2b<-ggplot(data = dts,aes(x=IOP))+geom_histogram()+facet_wrap(~adams_affected)
print(d1)## Warning: Removed 2 rows containing non-finite values (stat_density).
print(d2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing non-finite values (stat_bin).
print(d2b)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing non-finite values (stat_bin).
#check NA values
filter(dts,is.na(IOP))## # A tibble: 2 x 13
## Dog Gender Status DOB Date `Age in Days`
## <chr> <chr> <chr> <dttm> <dttm> <dbl>
## 1 Libe~ F achm ~ 2014-05-19 00:00:00 2014-07-22 00:00:00 64
## 2 Libe~ F achm ~ 2014-05-19 00:00:00 2014-07-22 00:00:00 64
## # ... with 7 more variables: `Age in Months` <dbl>, Time <chr>,
## # Tonomoter <chr>, Eye <chr>, IOP <dbl>, Puppy <chr>,
## # adams_affected <lgl>#confirmed missing data, delete
dts<-filter(dts,!is.na(IOP))
dim(dts)## [1] 6953 13#start descriptive analyses and graphics
#1: IOP trends over time, for three different times of measurement
p1<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_point()
print(p1)
p2<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_smooth()
print(p2)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p3<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_smooth()+facet_wrap(~Status)
print(p3)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p4<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_smooth()+facet_wrap(~adams_affected)
print(p4)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p4<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_smooth()+facet_wrap(~adams_affected)
print(p4)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p5<-ggplot(data = dts,aes(x=`Age in Days`,y=IOP,color=Time))+geom_smooth()+facet_wrap(~adams_affected+Gender)
print(p5)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
inferences by grouping data
dts<-mutate(dts,years=round(`Age in Days`/365))
dts## # A tibble: 6,953 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 6,943 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>summary(dts$years)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.826 3.000 7.000table(dts$years)##
## 0 1 2 3 4 5 6 7
## 1872 2007 1097 522 696 471 228 60Not Affected
#----------------------------
adams_dogs<-filter(dts,!adams_affected)
adams_dogs## # A tibble: 2,842 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 2,832 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>table(adams_dogs$years)##
## 0 1 2 3 4 5 6 7
## 838 684 444 288 204 144 180 60ftb<-table(adams_dogs$Dog,adams_dogs$years)
#number of dogs per age group
colSums(ftb>0)## 0 1 2 3 4 5 6 7
## 15 11 11 11 3 4 4 2# 0 1 2 3 4 5 6
#28 25 15 7 9 8 3
#raw mean by age group
adams_dogs_by_age<-group_by(adams_dogs,years)
adams_dogs_by_age## # A tibble: 2,842 x 14
## # Groups: years [8]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 2,832 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>summarize(adams_dogs_by_age,mn=mean(IOP),stdev=sd(IOP))## # A tibble: 8 x 3
## years mn stdev
## <dbl> <dbl> <dbl>
## 1 0 10.9 2.61
## 2 1 15.6 2.94
## 3 2 14.6 3.18
## 4 3 14.5 3.00
## 5 4 13.7 2.45
## 6 5 14.2 2.58
## 7 6 13.5 3.61
## 8 7 13.8 2.49adams_dogs_by_age<-group_by(adams_dogs,years,Dog)
adams_dogs_by_age## # A tibble: 2,842 x 14
## # Groups: years, Dog [61]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 2,832 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>mean_by_dog_age<-summarize(adams_dogs_by_age,mIOP=mean(IOP))
mean_by_dog_age<-ungroup(mean_by_dog_age)%>%group_by(years)
summarize(mean_by_dog_age,mn=mean(mIOP),stdev=sd(mIOP))## # A tibble: 8 x 3
## years mn stdev
## <dbl> <dbl> <dbl>
## 1 0 12.0 1.79
## 2 1 15.5 1.68
## 3 2 14.8 2.53
## 4 3 14.4 1.99
## 5 4 13.8 1.21
## 6 5 15.1 2.14
## 7 6 13.1 0.945
## 8 7 14.0 1.59g0<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_point()+geom_smooth()+ylim(5,30)
g1<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_smooth()+ylim(5,30)
print(g0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
print(g1)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#robust loess smoothing line and its SE BAND
ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_point()+geom_smooth(method = "lm")+ylim(5,30)
mean_by_dog_age<-mutate(mean_by_dog_age,cat_year=as.factor(years),
cyear=years-3,
cyear2=cyear^2,
cyear3=cyear^3)## Warning in mutate_impl(.data, dots, caller_env()): Unequal factor levels:
## coercing to character## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vectormean_by_dog_age## # A tibble: 61 x 7
## # Groups: years [8]
## years Dog mIOP cat_year cyear cyear2 cyear3
## <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 0 Arlington 9.49 0 -3 9 -27
## 2 0 Cleopatra 11.5 0 -3 9 -27
## 3 0 Esteban 12.7 0 -3 9 -27
## 4 0 Felix 12.7 0 -3 9 -27
## 5 0 Harry 13.8 0 -3 9 -27
## 6 0 Hope 10.2 0 -3 9 -27
## 7 0 Howgwarts 14.4 0 -3 9 -27
## 8 0 Javier 14.6 0 -3 9 -27
## 9 0 Jax 12.3 0 -3 9 -27
## 10 0 Liberty 10 0 -3 9 -27
## # ... with 51 more rowslm1<-lm(mIOP~Dog+cat_year,data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
aov(lm1)## Call:
## aov(formula = lm1)
##
## Terms:
## Dog cat_year Residuals
## Sum of Squares 210.34106 44.64018 45.16753
## Deg. of Freedom 30 7 23
##
## Residual standard error: 1.401358
## Estimated effects may be unbalancedanova(lm1)## Analysis of Variance Table
##
## Response: mIOP
## Df Sum Sq Mean Sq F value Pr(>F)
## Dog 30 210.341 7.0114 3.5703 0.001248 **
## cat_year 7 44.640 6.3772 3.2474 0.015118 *
## Residuals 23 45.168 1.9638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:tidyr':
##
## expandlm2<-lmer(mIOP~cat_year+(1|Dog),data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 232.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4407 -0.5851 -0.1109 0.5398 2.7365
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 1.185 1.088
## Residual 2.637 1.624
## Number of obs: 61, groups: Dog, 31
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.0888 0.4980 24.277
## cat_year1 3.3256 0.7313 4.547
## cat_year2 2.6811 0.7382 3.632
## cat_year3 2.2019 0.7391 2.979
## cat_year4 1.0399 1.1643 0.893
## cat_year5 2.3619 1.0480 2.254
## cat_year6 0.6852 1.0690 0.641
## cat_year7 1.2555 1.4180 0.885
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5 ct_yr6
## cat_year1 -0.624
## cat_year2 -0.632 0.508
## cat_year3 -0.632 0.496 0.551
## cat_year4 -0.417 0.286 0.342 0.369
## cat_year5 -0.466 0.314 0.366 0.389 0.355
## cat_year6 -0.461 0.299 0.328 0.333 0.260 0.326
## cat_year7 -0.349 0.223 0.235 0.238 0.174 0.256 0.333anova(lm2)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## cat_year 7 67.889 9.6984 3.6773library(lsmeans)## 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'.leastsquare<-lsmeans(lm1,"cat_year")
leastsquare2<-lsmeans(lm2,"cat_year")
leastsquare## cat_year lsmean SE df lower.CL upper.CL
## 0 13.51 0.635 23 12.20 14.8
## 1 15.24 0.628 23 13.94 16.5
## 2 14.03 0.631 23 12.72 15.3
## 3 13.41 0.636 23 12.09 14.7
## 4 10.80 1.106 23 8.51 13.1
## 5 11.69 1.044 23 9.53 13.9
## 6 9.15 1.327 23 6.40 11.9
## 7 8.81 1.684 23 5.33 12.3
##
## Results are averaged over the levels of: Dog
## Confidence level used: 0.95leastsquare2## cat_year lsmean SE df lower.CL upper.CL
## 0 12.1 0.503 50.9 11.1 13.1
## 1 15.4 0.583 53.0 14.2 16.6
## 2 14.8 0.583 53.0 13.6 15.9
## 3 14.3 0.584 53.0 13.1 15.5
## 4 13.1 1.087 46.7 10.9 15.3
## 5 14.5 0.952 49.9 12.5 16.4
## 6 12.8 0.966 53.0 10.8 14.7
## 7 13.3 1.364 52.2 10.6 16.1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(c0v1=c(-1,1,0,0,0,0,0,0),
c1v2=c(0,-1,1,0,0,0,0,0),
c2v3=c(0,0,-1,1,0,0,0,0),
c3v4=c(0,0,0,-1,1,0,0,0),
c4v5=c(0,0,0,0,-1,1,0,0),
c5v6=c(0,0,0,0,0,-1,1,0),
c6v7=c(0,0,0,0,0,-1,1,0))
ctrs## $c0v1
## [1] -1 1 0 0 0 0 0 0
##
## $c1v2
## [1] 0 -1 1 0 0 0 0 0
##
## $c2v3
## [1] 0 0 -1 1 0 0 0 0
##
## $c3v4
## [1] 0 0 0 -1 1 0 0 0
##
## $c4v5
## [1] 0 0 0 0 -1 1 0 0
##
## $c5v6
## [1] 0 0 0 0 0 -1 1 0
##
## $c6v7
## [1] 0 0 0 0 0 -1 1 0contrast(leastsquare,ctrs)## contrast estimate SE df t.ratio p.value
## c0v1 1.73 1.043 23 1.660 0.1106
## c1v2 -1.22 0.704 23 -1.726 0.0977
## c2v3 -0.62 0.612 23 -1.013 0.3214
## c3v4 -2.60 1.057 23 -2.465 0.0216
## c4v5 0.89 1.115 23 0.799 0.4326
## c5v6 -2.55 1.262 23 -2.017 0.0555
## c6v7 -2.55 1.262 23 -2.017 0.0555
##
## Results are averaged over the levels of: Dogdby<-table(mean_by_dog_age$Dog,mean_by_dog_age$years)
dby##
## 0 1 2 3 4 5 6 7
## Anjie 0 0 1 1 0 0 0 0
## Arlington 1 0 0 0 0 0 0 0
## Battle Creek 0 1 0 0 0 0 0 0
## Britney Spears 0 1 0 0 0 0 0 0
## Capone 0 0 1 1 1 1 0 0
## Cleopatra 1 1 0 0 0 0 0 0
## Esteban 1 0 0 0 0 0 0 0
## Felix 1 0 0 0 0 0 0 0
## Gilda 0 1 1 1 0 0 0 0
## Harry 1 1 1 1 0 0 0 0
## Hermione 0 0 0 0 0 0 1 0
## Hope 1 0 0 0 0 0 0 0
## Howgwarts 1 1 1 1 0 0 0 0
## Javier 1 0 0 0 0 0 0 0
## Jax 1 0 0 0 0 0 0 0
## Jayne 0 0 0 0 0 0 1 1
## Lacey 0 0 1 1 0 0 0 0
## Lady 0 1 1 1 0 0 0 0
## Liberty 1 0 0 0 0 0 0 0
## Macho 0 0 1 1 1 1 1 0
## Mateo 1 0 0 0 0 0 0 0
## Mattawan 0 1 1 0 0 0 0 0
## Nadia 0 0 0 0 0 1 1 1
## Piper 0 1 1 1 0 0 0 0
## Porkchop 0 1 0 0 0 0 0 0
## Red Poppy 1 0 0 0 0 0 0 0
## Spinning 0 0 0 1 1 1 0 0
## Stefano 1 0 0 0 0 0 0 0
## Turtle Dove 0 1 1 1 0 0 0 0
## Waterloo 1 0 0 0 0 0 0 0
## Yankee Doodle 1 0 0 0 0 0 0 0t(dby)%*%dby##
## 0 1 2 3 4 5 6 7
## 0 15 3 2 2 0 0 0 0
## 1 3 11 7 6 0 0 0 0
## 2 2 7 11 10 2 2 1 0
## 3 2 6 10 11 3 3 1 0
## 4 0 0 2 3 3 3 1 0
## 5 0 0 2 3 3 4 2 1
## 6 0 0 1 1 1 2 4 2
## 7 0 0 0 0 0 1 2 2nrow(dby) #number of dogs## [1] 31#31
table(mean_by_dog_age$Dog,mean_by_dog_age$years)%>%colSums()## 0 1 2 3 4 5 6 7
## 15 11 11 11 3 4 4 2leastsquare2## cat_year lsmean SE df lower.CL upper.CL
## 0 12.1 0.503 50.9 11.1 13.1
## 1 15.4 0.583 53.0 14.2 16.6
## 2 14.8 0.583 53.0 13.6 15.9
## 3 14.3 0.584 53.0 13.1 15.5
## 4 13.1 1.087 46.7 10.9 15.3
## 5 14.5 0.952 49.9 12.5 16.4
## 6 12.8 0.966 53.0 10.8 14.7
## 7 13.3 1.364 52.2 10.6 16.1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95g2<- ggplot(as.data.frame(leastsquare2), aes(x=cat_year, y=lsmean)) +
geom_bar(stat="identity", color="black",fill="darkgreen",
position=position_dodge()) +
geom_errorbar(aes(ymin=lsmean-SE, ymax=lsmean+SE), width=.2,
position=position_dodge(.9)) +ylim(0,35)
print(g2)
g3<- ggplot(as.data.frame(leastsquare2[-7,]), aes(x=cat_year, y=lsmean)) +
geom_bar(stat="identity", color="black",fill="darkgreen",
position=position_dodge()) +
geom_errorbar(aes(ymin=lsmean-SE, ymax=lsmean+SE), width=.2,
position=position_dodge(.9)) +ylim(0,35)
print(g3)
#compute overall year effect p-value.
#sex effect
#normal dogs
contrast(leastsquare2,ctrs)## contrast estimate SE df t.ratio p.value
## c0v1 3.326 0.745 51.6 4.465 <.0001
## c1v2 -0.645 0.739 35.3 -0.872 0.3889
## c2v3 -0.479 0.703 27.9 -0.682 0.5008
## c3v4 -1.162 1.148 34.9 -1.012 0.3184
## c4v5 1.322 1.267 28.6 1.044 0.3053
## c5v6 -1.677 1.251 38.5 -1.340 0.1880
## c6v7 -1.677 1.251 38.5 -1.340 0.1880summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 232.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4407 -0.5851 -0.1109 0.5398 2.7365
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 1.185 1.088
## Residual 2.637 1.624
## Number of obs: 61, groups: Dog, 31
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.0888 0.4980 24.277
## cat_year1 3.3256 0.7313 4.547
## cat_year2 2.6811 0.7382 3.632
## cat_year3 2.2019 0.7391 2.979
## cat_year4 1.0399 1.1643 0.893
## cat_year5 2.3619 1.0480 2.254
## cat_year6 0.6852 1.0690 0.641
## cat_year7 1.2555 1.4180 0.885
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5 ct_yr6
## cat_year1 -0.624
## cat_year2 -0.632 0.508
## cat_year3 -0.632 0.496 0.551
## cat_year4 -0.417 0.286 0.342 0.369
## cat_year5 -0.466 0.314 0.366 0.389 0.355
## cat_year6 -0.461 0.299 0.328 0.333 0.260 0.326
## cat_year7 -0.349 0.223 0.235 0.238 0.174 0.256 0.333#methods:
#pooled (averaged) data from 1200 to 1800 days of age and perfomed ANOVA by adams status
time_of_day<-filter(dts,`Age in Days`>=1200,`Age in Days`<=1800)%>%group_by(Dog,adams_affected,Time)%>%summarize(mIOP=mean(IOP))
time_of_day<-ungroup(time_of_day)
time_of_day## # A tibble: 36 x 4
## Dog adams_affected Time mIOP
## <chr> <lgl> <chr> <dbl>
## 1 Andy TRUE AM 27.6
## 2 Andy TRUE NOON 23.6
## 3 Andy TRUE PM 21.4
## 4 Angel TRUE AM 37.4
## 5 Angel TRUE NOON 29.3
## 6 Angel TRUE PM 26.4
## 7 Arnold TRUE AM 29.2
## 8 Arnold TRUE NOON 24.0
## 9 Arnold TRUE PM 22.1
## 10 Baron TRUE AM 22.8
## # ... with 26 more rowsexample<-filter(time_of_day,adams_affected==T)
library(lme4)
library(lsmeans)
lm1<-lmer(mIOP~Time+(1|Dog),data=example)
lscomp<-lsmeans(lm1,"Time")
lscomp## Time lsmean SE df lower.CL upper.CL
## AM 27.6 1.47 11.7 24.4 30.9
## NOON 24.1 1.47 11.7 20.8 27.3
## PM 23.0 1.47 11.7 19.8 26.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(NOON_AM=c(-1,1,0),
PM_NOON=c(0,-1,1),
PM_AM=c(-1,0,1))
contrast(lscomp,ctrs)## contrast estimate SE df t.ratio p.value
## NOON_AM -3.59 1.09 16 -3.308 0.0044
## PM_NOON -1.06 1.09 16 -0.978 0.3424
## PM_AM -4.66 1.09 16 -4.286 0.0006example<-filter(time_of_day,adams_affected==F)
library(lme4)
library(lsmeans)
lm1<-lmer(mIOP~Time+(1|Dog),data=example)
lscomp<-lsmeans(lm1,"Time")
lscomp## Time lsmean SE df lower.CL upper.CL
## AM 14.5 0.922 2.09 10.69 18.3
## NOON 14.0 0.922 2.09 10.16 17.8
## PM 13.5 0.922 2.09 9.73 17.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(NOON_AM=c(-1,1,0),
PM_NOON=c(0,-1,1),
PM_AM=c(-1,0,1))
contrast(lscomp,ctrs)## contrast estimate SE df t.ratio p.value
## NOON_AM -0.530 0.238 4 -2.232 0.0895
## PM_NOON -0.433 0.238 4 -1.823 0.1424
## PM_AM -0.963 0.238 4 -4.054 0.0154Not Affected, dogs younger than 1 y.o. are deleted
#----------------------------
adams_dogs<-filter(dts,!adams_affected,`Age in Days`>=366)
adams_dogs## # A tibble: 1,596 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 1,586 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>table(adams_dogs$years)##
## 1 2 3 4 5 6 7
## 276 444 288 204 144 180 60ftb<-table(adams_dogs$Dog,adams_dogs$years)
#number of dogs per age group
colSums(ftb>0)## 1 2 3 4 5 6 7
## 9 11 11 3 4 4 2print(ftb)##
## 1 2 3 4 5 6 7
## Anjie 0 12 6 0 0 0 0
## Battle Creek 24 0 0 0 0 0 0
## Capone 0 36 42 66 60 0 0
## Gilda 36 60 30 0 0 0 0
## Harry 36 36 6 0 0 0 0
## Hermione 0 0 0 0 0 18 0
## Howgwarts 42 36 6 0 0 0 0
## Jayne 0 0 0 0 0 78 24
## Lacey 0 30 6 0 0 0 0
## Lady 36 72 30 0 0 0 0
## Macho 0 6 66 72 42 12 0
## Mattawan 6 24 0 0 0 0 0
## Nadia 0 0 0 0 6 72 36
## Piper 30 66 36 0 0 0 0
## Porkchop 36 0 0 0 0 0 0
## Spinning 0 0 24 66 36 0 0
## Turtle Dove 30 66 36 0 0 0 0#raw mean by age group
adams_dogs_by_age<-group_by(adams_dogs,years)
adams_dogs_by_age## # A tibble: 1,596 x 14
## # Groups: years [7]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 1,586 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>adams_dogs_by_age<-group_by(adams_dogs,years,Dog)
mean_by_dog_age<-summarize(adams_dogs_by_age,mIOP=mean(IOP))
mean_by_dog_age<-ungroup(mean_by_dog_age)%>%group_by(years)
summarize(mean_by_dog_age,mn=mean(mIOP),stdev=sd(mIOP))## # A tibble: 7 x 3
## years mn stdev
## <dbl> <dbl> <dbl>
## 1 1 16.8 2.77
## 2 2 14.8 2.53
## 3 3 14.4 1.99
## 4 4 13.8 1.21
## 5 5 15.1 2.14
## 6 6 13.1 0.945
## 7 7 14.0 1.59g0<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_point()+geom_smooth()+ylim(5,30)
g1<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_smooth()+ylim(5,30)
print(g0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
print(g1)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Affected
#----------------------------
adams_dogs<-filter(dts,adams_affected)
adams_dogs## # A tibble: 4,111 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 4,101 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>table(adams_dogs$years)##
## 0 1 2 3 4 5 6
## 1034 1323 653 234 492 327 48ftb<-table(adams_dogs$Dog,adams_dogs$years)
#number of dogs per age group
colSums(ftb>0)## 0 1 2 3 4 5 6
## 28 25 15 7 9 8 3# 0 1 2 3 4 5 6
#28 25 15 7 9 8 3
#raw mean by age group
adams_dogs_by_age<-group_by(adams_dogs,years)
adams_dogs_by_age## # A tibble: 4,111 x 14
## # Groups: years [7]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 4,101 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>summarize(adams_dogs_by_age,mn=mean(IOP),stdev=sd(IOP))## # A tibble: 7 x 3
## years mn stdev
## <dbl> <dbl> <dbl>
## 1 0 12.8 3.55
## 2 1 18.0 3.31
## 3 2 19.9 4.30
## 4 3 23.1 6.25
## 5 4 24.7 6.15
## 6 5 24.2 7.08
## 7 6 28.1 5.73adams_dogs_by_age<-group_by(adams_dogs,years,Dog)
adams_dogs_by_age## # A tibble: 4,111 x 14
## # Groups: years, Dog [95]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 4,101 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>mean_by_dog_age<-summarize(adams_dogs_by_age,mIOP=mean(IOP))
mean_by_dog_age<-ungroup(mean_by_dog_age)%>%group_by(years)
summarize(mean_by_dog_age,mn=mean(mIOP),stdev=sd(mIOP))## # A tibble: 7 x 3
## years mn stdev
## <dbl> <dbl> <dbl>
## 1 0 13.8 2.19
## 2 1 18.5 1.99
## 3 2 20.5 3.44
## 4 3 24.1 4.23
## 5 4 24.7 3.63
## 6 5 25.0 4.47
## 7 6 27.9 3.41g0<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_point()+geom_smooth()+ylim(5,30)
g1<-ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_smooth()+ylim(5,30)
print(g0)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 4 rows containing non-finite values (stat_smooth).## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.03## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.03## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.2668e-016## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.03## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.03## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 1.2668e-016## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 1## Warning: Removed 4 rows containing missing values (geom_point).
print(g1)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 4 rows containing non-finite values (stat_smooth).## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.03## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.03## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.2668e-016## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.03## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.03## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 1.2668e-016## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 1
#robust loess smoothing line and its SE BAND
ggplot(data = mean_by_dog_age,
aes(x=`years`,y=mIOP))+geom_point()+geom_smooth(method = "lm")+ylim(5,30)## Warning: Removed 4 rows containing non-finite values (stat_smooth).## Warning: Removed 4 rows containing missing values (geom_point).
mean_by_dog_age<-mutate(mean_by_dog_age,cat_year=as.factor(years),
cyear=years-3,
cyear2=cyear^2,
cyear3=cyear^3)## Warning in mutate_impl(.data, dots, caller_env()): Unequal factor levels:
## coercing to character## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vectormean_by_dog_age## # A tibble: 95 x 7
## # Groups: years [7]
## years Dog mIOP cat_year cyear cyear2 cyear3
## <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 0 Aldo 17.6 0 -3 9 -27
## 2 0 Anja 17.6 0 -3 9 -27
## 3 0 Anton 15.8 0 -3 9 -27
## 4 0 Armin 16.9 0 -3 9 -27
## 5 0 Betty 14.6 0 -3 9 -27
## 6 0 Czar 14.6 0 -3 9 -27
## 7 0 Dandelion 10.4 0 -3 9 -27
## 8 0 Daphne 11.5 0 -3 9 -27
## 9 0 Dogwood 11.8 0 -3 9 -27
## 10 0 Dryas 11.4 0 -3 9 -27
## # ... with 85 more rowslm1<-lm(mIOP~Dog+cat_year,data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
aov(lm1)## Call:
## aov(formula = lm1)
##
## Terms:
## Dog cat_year Residuals
## Sum of Squares 1701.6606 649.6354 210.8523
## Deg. of Freedom 36 6 52
##
## Residual standard error: 2.013666
## Estimated effects may be unbalancedanova(lm1)## Analysis of Variance Table
##
## Response: mIOP
## Df Sum Sq Mean Sq F value Pr(>F)
## Dog 36 1701.66 47.268 11.657 7.92e-15 ***
## cat_year 6 649.64 108.273 26.702 2.91e-14 ***
## Residuals 52 210.85 4.055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lme4)
lm2<-lmer(mIOP~cat_year+(1|Dog),data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 440.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9462 -0.4696 -0.1499 0.3180 2.3099
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 4.750 2.18
## Residual 4.325 2.08
## Number of obs: 95, groups: Dog, 37
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.1773 0.5560 25.499
## cat_year1 4.8068 0.5807 8.278
## cat_year2 7.0293 0.7083 9.924
## cat_year3 9.2587 1.0748 8.614
## cat_year4 9.8631 1.0338 9.541
## cat_year5 11.0264 1.0482 10.520
## cat_year6 15.2756 1.4873 10.271
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5
## cat_year1 -0.481
## cat_year2 -0.418 0.398
## cat_year3 -0.403 0.255 0.271
## cat_year4 -0.449 0.264 0.253 0.442
## cat_year5 -0.424 0.261 0.253 0.389 0.504
## cat_year6 -0.289 0.185 0.179 0.241 0.327 0.345anova(lm2)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## cat_year 6 974.55 162.42 37.55library(lsmeans)
leastsquare<-lsmeans(lm1,"cat_year")
leastsquare2<-lsmeans(lm2,"cat_year")
leastsquare## cat_year lsmean SE df lower.CL upper.CL
## 0 14.8 0.471 52 13.9 15.8
## 1 19.7 0.498 52 18.7 20.7
## 2 21.8 0.613 52 20.5 23.0
## 3 21.9 1.061 52 19.7 24.0
## 4 22.2 1.092 52 20.1 24.4
## 5 23.9 1.074 52 21.8 26.1
## 6 28.9 1.492 52 25.9 31.9
##
## Results are averaged over the levels of: Dog
## Confidence level used: 0.95leastsquare2## cat_year lsmean SE df lower.CL upper.CL
## 0 14.2 0.558 68.9 13.1 15.3
## 1 19.0 0.582 73.5 17.8 20.1
## 2 21.2 0.698 87.3 19.8 22.6
## 3 23.4 1.001 87.4 21.4 25.4
## 4 24.0 0.939 85.6 22.2 25.9
## 5 25.2 0.966 87.8 23.3 27.1
## 6 29.5 1.442 77.5 26.6 32.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(c0v1=c(-1,1,0,0,0,0,0),
c1v2=c(0,-1,1,0,0,0,0),
c2v3=c(0,0,-1,1,0,0,0),
c3v4=c(0,0,0,-1,1,0,0),
c4v5=c(0,0,0,0,-1,1,0),
c5v6=c(0,0,0,0,0,-1,1))
ctrs## $c0v1
## [1] -1 1 0 0 0 0 0
##
## $c1v2
## [1] 0 -1 1 0 0 0 0
##
## $c2v3
## [1] 0 0 -1 1 0 0 0
##
## $c3v4
## [1] 0 0 0 -1 1 0 0
##
## $c4v5
## [1] 0 0 0 0 -1 1 0
##
## $c5v6
## [1] 0 0 0 0 0 -1 1contrast(leastsquare,ctrs)## contrast estimate SE df t.ratio p.value
## c0v1 4.8401 0.570 52 8.498 <.0001
## c1v2 2.1099 0.712 52 2.962 0.0046
## c2v3 0.0868 1.287 52 0.067 0.9465
## c3v4 0.3830 1.121 52 0.342 0.7341
## c4v5 1.6828 1.018 52 1.653 0.1044
## c5v6 5.0165 1.479 52 3.392 0.0013
##
## Results are averaged over the levels of: Dogdby<-table(mean_by_dog_age$Dog,mean_by_dog_age$years)
dby##
## 0 1 2 3 4 5 6
## Aldo 1 1 1 0 0 0 0
## Andy 0 0 0 1 1 0 0
## Angel 0 0 0 0 1 1 0
## Anja 1 1 1 0 0 0 0
## Anton 1 1 1 0 0 0 0
## Armin 1 1 1 1 0 1 1
## Arnold 0 0 0 0 1 1 0
## Baron 0 0 0 0 1 1 1
## Betty 1 1 0 0 0 0 0
## Chubbie 0 0 0 0 1 1 1
## Czar 1 1 1 1 1 1 0
## Dandelion 1 1 1 0 0 0 0
## Daphne 1 1 0 0 0 0 0
## Dogwood 1 1 0 0 0 0 0
## Dryas 1 1 0 0 0 0 0
## Emulsion 1 1 1 0 0 0 0
## Fargo 1 1 0 0 0 0 0
## Gaston 1 1 1 0 0 0 0
## Gepetto 1 1 1 0 0 0 0
## Goldilocks 1 1 1 0 0 0 0
## Granny 1 1 1 0 0 0 0
## Gretel 1 1 1 0 0 0 0
## Grimm 1 1 1 0 0 0 0
## Grumpy 1 1 1 0 0 0 0
## Hoss 0 0 1 1 0 0 0
## Kepler 1 1 0 0 0 0 0
## Kuiper 1 0 0 0 0 0 0
## Laser 1 1 0 0 0 0 0
## Lathe 1 1 0 0 0 0 0
## Lepton 1 1 0 0 0 0 0
## Lumi 1 1 0 0 0 0 0
## Luminous 1 0 0 0 0 0 0
## Lunar 1 0 0 0 0 0 0
## Mabel 1 1 0 0 0 0 0
## Princess 0 0 0 1 1 0 0
## Rita 0 0 0 1 1 1 0
## Toby 0 0 0 1 1 1 0t(dby)%*%dby##
## 0 1 2 3 4 5 6
## 0 28 25 14 2 1 2 1
## 1 25 25 14 2 1 2 1
## 2 14 14 15 3 1 2 1
## 3 2 2 3 7 5 4 1
## 4 1 1 1 5 9 7 2
## 5 2 2 2 4 7 8 3
## 6 1 1 1 1 2 3 3nrow(dby) #number of dogs## [1] 37#31
table(mean_by_dog_age$Dog,mean_by_dog_age$years)%>%colSums()## 0 1 2 3 4 5 6
## 28 25 15 7 9 8 3leastsquare2## cat_year lsmean SE df lower.CL upper.CL
## 0 14.2 0.558 68.9 13.1 15.3
## 1 19.0 0.582 73.5 17.8 20.1
## 2 21.2 0.698 87.3 19.8 22.6
## 3 23.4 1.001 87.4 21.4 25.4
## 4 24.0 0.939 85.6 22.2 25.9
## 5 25.2 0.966 87.8 23.3 27.1
## 6 29.5 1.442 77.5 26.6 32.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95g2<- ggplot(as.data.frame(leastsquare2), aes(x=cat_year, y=lsmean)) +
geom_bar(stat="identity", color="black",fill="darkgreen",
position=position_dodge()) +
geom_errorbar(aes(ymin=lsmean-SE, ymax=lsmean+SE), width=.2,
position=position_dodge(.9))+ylim(0,35)
print(g2)
g3<- ggplot(as.data.frame(leastsquare2[-7,]), aes(x=cat_year, y=lsmean)) +
geom_bar(stat="identity", color="black",fill="darkgreen",
position=position_dodge()) +
geom_errorbar(aes(ymin=lsmean-SE, ymax=lsmean+SE), width=.2,
position=position_dodge(.9)) +ylim(0,35)
print(g3)
#compute overall year effect p-value.
#sex effect
#normal dogs
contrast(leastsquare2,ctrs)## contrast estimate SE df t.ratio p.value
## c0v1 4.807 0.581 54.8 8.268 <.0001
## c1v2 2.222 0.718 57.9 3.095 0.0030
## c2v3 2.229 1.127 76.2 1.977 0.0516
## c3v4 0.604 1.119 59.2 0.540 0.5912
## c4v5 1.163 1.039 55.2 1.119 0.2678
## c5v6 4.249 1.501 57.4 2.831 0.0064summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 440.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9462 -0.4696 -0.1499 0.3180 2.3099
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 4.750 2.18
## Residual 4.325 2.08
## Number of obs: 95, groups: Dog, 37
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.1773 0.5560 25.499
## cat_year1 4.8068 0.5807 8.278
## cat_year2 7.0293 0.7083 9.924
## cat_year3 9.2587 1.0748 8.614
## cat_year4 9.8631 1.0338 9.541
## cat_year5 11.0264 1.0482 10.520
## cat_year6 15.2756 1.4873 10.271
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5
## cat_year1 -0.481
## cat_year2 -0.418 0.398
## cat_year3 -0.403 0.255 0.271
## cat_year4 -0.449 0.264 0.253 0.442
## cat_year5 -0.424 0.261 0.253 0.389 0.504
## cat_year6 -0.289 0.185 0.179 0.241 0.327 0.345#methods:
#pooled (averaged) data from 1200 to 1800 days of age and perfomed ANOVA by adams status
time_of_day<-filter(dts,`Age in Days`>=1200,`Age in Days`<=1800)%>%group_by(Dog,adams_affected,Time)%>%summarize(mIOP=mean(IOP))
time_of_day<-ungroup(time_of_day)
time_of_day## # A tibble: 36 x 4
## Dog adams_affected Time mIOP
## <chr> <lgl> <chr> <dbl>
## 1 Andy TRUE AM 27.6
## 2 Andy TRUE NOON 23.6
## 3 Andy TRUE PM 21.4
## 4 Angel TRUE AM 37.4
## 5 Angel TRUE NOON 29.3
## 6 Angel TRUE PM 26.4
## 7 Arnold TRUE AM 29.2
## 8 Arnold TRUE NOON 24.0
## 9 Arnold TRUE PM 22.1
## 10 Baron TRUE AM 22.8
## # ... with 26 more rowsexample<-filter(time_of_day,adams_affected==T)
library(lme4)
library(lsmeans)
lm1<-lmer(mIOP~Time+(1|Dog),data=example)
lscomp<-lsmeans(lm1,"Time")
lscomp## Time lsmean SE df lower.CL upper.CL
## AM 27.6 1.47 11.7 24.4 30.9
## NOON 24.1 1.47 11.7 20.8 27.3
## PM 23.0 1.47 11.7 19.8 26.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(NOON_AM=c(-1,1,0),
PM_NOON=c(0,-1,1),
PM_AM=c(-1,0,1))
contrast(lscomp,ctrs)## contrast estimate SE df t.ratio p.value
## NOON_AM -3.59 1.09 16 -3.308 0.0044
## PM_NOON -1.06 1.09 16 -0.978 0.3424
## PM_AM -4.66 1.09 16 -4.286 0.0006example<-filter(time_of_day,adams_affected==F)
library(lme4)
library(lsmeans)
lm1<-lmer(mIOP~Time+(1|Dog),data=example)
lscomp<-lsmeans(lm1,"Time")
lscomp## Time lsmean SE df lower.CL upper.CL
## AM 14.5 0.922 2.09 10.69 18.3
## NOON 14.0 0.922 2.09 10.16 17.8
## PM 13.5 0.922 2.09 9.73 17.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(NOON_AM=c(-1,1,0),
PM_NOON=c(0,-1,1),
PM_AM=c(-1,0,1))
contrast(lscomp,ctrs)## contrast estimate SE df t.ratio p.value
## NOON_AM -0.530 0.238 4 -2.232 0.0895
## PM_NOON -0.433 0.238 4 -1.823 0.1424
## PM_AM -0.963 0.238 4 -4.054 0.0154gender effect affected
#----------------------------
adams_dogs<-filter(dts,adams_affected)
adams_dogs## # A tibble: 4,111 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 4,101 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>ftb<-table(adams_dogs$Dog,adams_dogs$years)
#number of dogs per age group
colSums(ftb>0)## 0 1 2 3 4 5 6
## 28 25 15 7 9 8 3# 0 1 2 3 4 5 6
#28 25 15 7 9 8 3
#raw mean by age group
adams_dogs_by_age<-group_by(adams_dogs,years,Dog,Gender)
adams_dogs_by_age## # A tibble: 4,111 x 14
## # Groups: years, Dog, Gender [95]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 2 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 3 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 4 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 5 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 6 Aldo M adamt~ 2012-03-02 00:00:00 2012-05-29 00:00:00
## 7 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 8 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 9 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## 10 Aldo M adamt~ 2012-03-02 00:00:00 2012-06-19 00:00:00
## # ... with 4,101 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>mean_by_dog_age<-summarize(adams_dogs_by_age,mIOP=mean(IOP))
mean_by_dog_age<-ungroup(mean_by_dog_age)%>%group_by(years)
mean_by_dog_age## # A tibble: 95 x 4
## # Groups: years [7]
## years Dog Gender mIOP
## <dbl> <chr> <chr> <dbl>
## 1 0 Aldo M 17.6
## 2 0 Anja F 17.6
## 3 0 Anton M 15.8
## 4 0 Armin M 16.9
## 5 0 Betty F 14.6
## 6 0 Czar M 14.6
## 7 0 Dandelion F 10.4
## 8 0 Daphne F 11.5
## 9 0 Dogwood M 11.8
## 10 0 Dryas M 11.4
## # ... with 85 more rowslibrary(lme4)
mean_by_dog_age<-mutate(mean_by_dog_age,cat_year=as.factor(years),
cyear=years-3,
cyear2=cyear^2,
cyear3=cyear^3)## Warning in mutate_impl(.data, dots, caller_env()): Unequal factor levels:
## coercing to character## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vectorlm2<-lmer(mIOP~cat_year+Gender+(1|Dog),data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + Gender + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 435.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8759 -0.5444 -0.1231 0.4082 2.5059
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 4.139 2.034
## Residual 4.412 2.100
## Number of obs: 95, groups: Dog, 37
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.9001 0.6836 21.797
## cat_year1 4.8488 0.5863 8.270
## cat_year2 7.1321 0.7154 9.969
## cat_year3 9.5720 1.0745 8.908
## cat_year4 10.1799 1.0278 9.904
## cat_year5 11.3225 1.0466 10.819
## cat_year6 15.6082 1.4978 10.421
## GenderM -1.4973 0.8250 -1.815
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5 ct_yr6
## cat_year1 -0.368
## cat_year2 -0.289 0.400
## cat_year3 -0.260 0.261 0.276
## cat_year4 -0.293 0.272 0.261 0.432
## cat_year5 -0.266 0.268 0.261 0.382 0.491
## cat_year6 -0.154 0.190 0.187 0.240 0.322 0.340
## GenderM -0.610 -0.044 -0.085 -0.101 -0.108 -0.119 -0.123anova(lm2)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## cat_year 6 1013.06 168.844 38.2719
## Gender 1 14.53 14.531 3.2937library(lsmeans)
leastsquare2<-lsmeans(lm2,"Gender")
leastsquare2## Gender lsmean SE df lower.CL upper.CL
## F 23.3 0.701 45.8 21.9 24.7
## M 21.8 0.579 35.9 20.6 23.0
##
## Results are averaged over the levels of: cat_year
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(gender=c(-1,1))
contrast(leastsquare2,ctrs)## contrast estimate SE df t.ratio p.value
## gender -1.5 0.826 34.9 -1.812 0.0786
##
## Results are averaged over the levels of: cat_yeargender effect not affected
#----------------------------
adams_dogs<-filter(dts,!adams_affected)
adams_dogs## # A tibble: 2,842 x 14
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 2,832 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>ftb<-table(adams_dogs$Dog,adams_dogs$years)
#number of dogs per age group
colSums(ftb>0)## 0 1 2 3 4 5 6 7
## 15 11 11 11 3 4 4 2# 0 1 2 3 4 5 6
#28 25 15 7 9 8 3
#raw mean by age group
adams_dogs_by_age<-group_by(adams_dogs,years,Dog,Gender)
adams_dogs_by_age## # A tibble: 2,842 x 14
## # Groups: years, Dog, Gender [61]
## Dog Gender Status DOB Date
## <chr> <chr> <chr> <dttm> <dttm>
## 1 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 2 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 3 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 4 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 5 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 6 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-06 00:00:00
## 7 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 8 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 9 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## 10 Anjie F achm ~ 2013-03-27 00:00:00 2015-03-12 00:00:00
## # ... with 2,832 more rows, and 9 more variables: `Age in Days` <dbl>,
## # `Age in Months` <dbl>, Time <chr>, Tonomoter <chr>, Eye <chr>,
## # IOP <dbl>, Puppy <chr>, adams_affected <lgl>, years <dbl>mean_by_dog_age<-summarize(adams_dogs_by_age,mIOP=mean(IOP))
mean_by_dog_age<-ungroup(mean_by_dog_age)%>%group_by(years)
mean_by_dog_age## # A tibble: 61 x 4
## # Groups: years [8]
## years Dog Gender mIOP
## <dbl> <chr> <chr> <dbl>
## 1 0 Arlington M 9.49
## 2 0 Cleopatra F 11.5
## 3 0 Esteban M 12.7
## 4 0 Felix M 12.7
## 5 0 Harry M 13.8
## 6 0 Hope F 10.2
## 7 0 Howgwarts F 14.4
## 8 0 Javier M 14.6
## 9 0 Jax M 12.3
## 10 0 Liberty F 10
## # ... with 51 more rowslibrary(lme4)
mean_by_dog_age<-mutate(mean_by_dog_age,cat_year=as.factor(years),
cyear=years-3,
cyear2=cyear^2,
cyear3=cyear^3)## Warning in mutate_impl(.data, dots, caller_env()): Unequal factor levels:
## coercing to character## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vector
## Warning in mutate_impl(.data, dots, caller_env()): binding character and
## factor vector, coercing into character vectorlm2<-lmer(mIOP~cat_year+Gender+(1|Dog),data=mean_by_dog_age)
#linear mixed model with randome effect of dog and a year-specific mean
summary(lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: mIOP ~ cat_year + Gender + (1 | Dog)
## Data: mean_by_dog_age
##
## REML criterion at convergence: 231.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45740 -0.53212 -0.05965 0.51106 2.64293
##
## Random effects:
## Groups Name Variance Std.Dev.
## Dog (Intercept) 1.151 1.073
## Residual 2.682 1.638
## Number of obs: 61, groups: Dog, 31
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.7827 0.6229 18.916
## cat_year1 3.4225 0.7431 4.606
## cat_year2 2.7960 0.7529 3.714
## cat_year3 2.3169 0.7536 3.074
## cat_year4 1.0180 1.1715 0.869
## cat_year5 2.4056 1.0528 2.285
## cat_year6 0.8857 1.0954 0.809
## cat_year7 1.5220 1.4530 1.048
## GenderM 0.5081 0.6287 0.808
##
## Correlation of Fixed Effects:
## (Intr) ct_yr1 ct_yr2 ct_yr3 ct_yr4 ct_yr5 ct_yr6 ct_yr7
## cat_year1 -0.586
## cat_year2 -0.605 0.519
## cat_year3 -0.604 0.507 0.561
## cat_year4 -0.305 0.274 0.325 0.351
## cat_year5 -0.386 0.313 0.362 0.384 0.349
## cat_year6 -0.484 0.320 0.351 0.355 0.242 0.320
## cat_year7 -0.394 0.246 0.261 0.264 0.159 0.252 0.355
## GenderM -0.598 0.153 0.177 0.175 -0.048 0.022 0.204 0.201anova(lm2)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## cat_year 7 68.701 9.8144 3.6592
## Gender 1 1.752 1.7520 0.6532library(lsmeans)
leastsquare2<-lsmeans(lm2,"Gender")
leastsquare2## Gender lsmean SE df lower.CL upper.CL
## F 13.6 0.465 27.9 12.6 14.5
## M 14.1 0.534 22.8 13.0 15.2
##
## Results are averaged over the levels of: cat_year
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95ctrs<-list(gender=c(-1,1))
contrast(leastsquare2,ctrs)## contrast estimate SE df t.ratio p.value
## gender 0.508 0.637 26.6 0.797 0.4323
##
## Results are averaged over the levels of: cat_year