Figures for “Assessment of Menstrual Health Status and Evolution through Mobile Apps for Fertility Awareness”"
Laura Symul, Katarzyna Wac, Paula Hillard and Marcel Salathé
2019
source("Scripts/libraries.R")##
## Attaching package: 'ggridges'## The following object is masked from 'package:ggplot2':
##
## scale_discrete_manual##
## Attaching package: 'lubridate'## The following objects are masked from 'package:data.table':
##
## hour, isoweek, mday, minute, month, quarter, second, wday,
## week, yday, year## The following object is masked from 'package:base':
##
## date##
## Attaching package: 'dplyr'## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union##
## Attaching package: 'reshape'## The following object is masked from 'package:dplyr':
##
## rename## The following object is masked from 'package:lubridate':
##
## stamp## The following object is masked from 'package:data.table':
##
## melt## Loading required package: foreach## Loading required package: iterators##
## Attaching package: 'plotrix'## The following object is masked from 'package:scales':
##
## rescale##
## Attaching package: 'chron'## The following object is masked from 'package:foreach':
##
## times## The following objects are masked from 'package:lubridate':
##
## days, hours, minutes, seconds, yearssource("Scripts/variables.R")
source("Scripts/functions.R")1 User history examples
1.1 History Kindara user
source("Scripts/_setup_Kindara.R")load(paste0(IO$kindara_data_days,'02 processed/day104-110.Rdata'), verbose = TRUE)## Loading objects:
## daysuser_id = '0a67d988dae24b669ab0644a656191d7'
j = (days$user_id == user_id)
start.date = as.Date('2015-01-01')
j = which(j & (days$date >= start.date))
d = days[j,]
plot.tracking.history(d = d)
pdf(paste0(IO$panels,'1A_history_kindara_example_',user_id,'.pdf'), useDingbats = FALSE, width = viz$pdf.full.width*0.75, height = viz$pdf.full.width/5)
plot.tracking.history(d = d)
dev.off()## quartz_off_screen
## 21.2 History Sympto user
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$restricted_figure_data,"history_sympto.Rdata"), verbose = TRUE)## Loading objects:
## cycledays_sympto_user_historyplot.tracking.history(d = cycledays_sympto_user_history)
u = unique(cycledays_sympto_user_history$user_id)
pdf(paste0(IO$panels,'1B_history_sympto_example_',u,'.pdf'), useDingbats = FALSE, width = viz$pdf.full.width*0.75, height = viz$pdf.full.width/5)
plot.tracking.history(d = cycledays_sympto_user_history, show_goal = TRUE, relative_date = TRUE)
dev.off()## quartz_off_screen
## 2rm(cycledays_sympto_user_history, u)2 Demographics (Fig 2A, S1C)
#sympto
load(file = paste0(IO$restricted_figure_data,"users_demo_sympto.Rdata"), verbose = TRUE)## Loading objects:
## users_demo_sympto# kindara
load(file = paste0(IO$restricted_figure_data,"users_demo_kindara.Rdata"), verbose = TRUE)## Loading objects:
## users_demo_kindara2.1 Sympto
summary(users_demo_sympto)## age_registration age_now menarche_year cm
## Min. : 0.00 Min. : 0.00 Min. : 8.00 Min. :100.0
## 1st Qu.:26.00 1st Qu.:29.00 1st Qu.:12.00 1st Qu.:160.0
## Median :29.00 Median :33.00 Median :13.00 Median :165.0
## Mean :29.94 Mean :33.52 Mean :12.85 Mean :165.3
## 3rd Qu.:34.00 3rd Qu.:38.00 3rd Qu.:14.00 3rd Qu.:170.0
## Max. :74.00 Max. :77.00 Max. :18.00 Max. :229.0
## NA's :2570 NA's :2570 NA's :2819 NA's :2735
## kg bmi cm_o bmi_o
## Min. : 20.00 Min. : 7.703 Min. : 50.0 Min. : 7.703
## 1st Qu.: 54.00 1st Qu.: 19.818 1st Qu.:160.0 1st Qu.:19.818
## Median : 60.00 Median : 21.631 Median :165.0 Median :21.630
## Mean : 62.23 Mean : 22.788 Mean :161.3 Mean :22.840
## 3rd Qu.: 67.00 3rd Qu.: 24.342 3rd Qu.:170.0 3rd Qu.:24.314
## Max. :149.00 Max. :103.939 Max. :229.0 Max. :86.565
## NA's :2765 NA's :2786 NA's :2735 NA's :3214round(apply(users_demo_sympto, 2, sd, na.rm = TRUE), digits = 2)## age_registration age_now menarche_year cm
## 6.36 6.69 1.61 6.73
## kg bmi cm_o bmi_o
## 13.14 4.76 21.25 5.02100 - round(apply(users_demo_sympto, 2, function(x) sum(is.na(x)))/nrow(users_demo_sympto)*100)## age_registration age_now menarche_year cm
## 81 81 79 80
## kg bmi cm_o bmi_o
## 80 80 80 762.2 Kindara
round(mean(users_demo_kindara$age_registration, na.rm = TRUE), digits = 1)## [1] 29.2round(sd(users_demo_kindara$age_registration, na.rm = TRUE), digits = 1)## [1] 5.9round(100 - (sum(is.na(users_demo_kindara$age_registration))/nrow(users_demo_kindara)*100))## [1] 132.3 Figures
plot_demographics_histogram = function(x, xlim = c(10, 60), ylim = NA, xlab = "x"){
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
breaks = min(x, na.rm = TRUE):(max(x, na.rm = TRUE)+1)-0.5
breaks = round(min(x, na.rm = TRUE)):round((max(x, na.rm = TRUE)+1))-0.5
h = hist(x, breaks = breaks, plot = FALSE)
if(any(is.na(ylim))){ylim = c(0,max(h$density)*1.03)}
plot(h$mids,h$counts/sum(h$counts), type = 'n', axes = FALSE, xlim = xlim, ylim = ylim,yaxs = 'i', xaxs = 'i',
xlab = xlab, ylab = '% of users', main = '')
at.x = seq(0,xlim[2]*1.5,by = 10)
axis(1, col = 'lightgray', col.ticks = 'lightgray')
at.y = seq(0,max(h$density),by = 0.01)
while(length(at.y) >= 10){at.y = at.y[rep(c(TRUE,FALSE))]}
axis(2, at = at.y , labels = 100*at.y, col = 'lightgray', col.ticks = 'lightgray')
abline(h = at.y, col = 'lightgray')
abline(v = at.x, col = 'lightgray')
polygon(c(h$mids,rev(h$mids)), c(h$counts/sum(h$counts),h$counts*0), col = viz_cols$colors_lsy[1], border = NA)
abline(v = median(x, na.rm = TRUE), lty = 2)
h = data.frame(x = h$mids, y = h$counts/sum(h$counts))
return(h)
}x2 = users_demo_kindara$age_now
h2 = plot_demographics_histogram(x = x2, xlim = c(10, 60), xlab = "Age at time of study")
x = users_demo_sympto$age_now
h = plot_demographics_histogram(x = x, xlim = c(10,60), ylim = c(0,0.08),xlab = 'Age at time of study')
points(h2$x, h2$y, type = "l", col = "mediumseagreen")
pdf.file = paste0(IO$panels,'2_demographics_age_now.pdf')
pdf(pdf.file, width = 3.5, height = 2)
h = plot_demographics_histogram(x = x, xlim = c(10,60), ylim = c(0,0.08),xlab = 'Age at time of study')
points(h2$x, h2$y, type = "l", col = "mediumseagreen")
dev.off()## quartz_off_screen
## 2x2 = users_demo_kindara$age_registration
h2 = plot_demographics_histogram(x = x2, xlim = c(10, 60), xlab = "Age at registration")
x = users_demo_sympto$age_registration
h = plot_demographics_histogram(x = x, xlim = c(10,60), ylim = c(0,0.08) ,xlab = 'Age at registration')
points(h2$x, h2$y, type = "l", col = "mediumseagreen")
pdf.file = paste0(IO$panels,'2_demographics_age_registration.pdf')
pdf(pdf.file, width = 3.5, height = 2)
h = plot_demographics_histogram(x = x, xlim = c(10,60), ylim = c(0,0.08),xlab = 'Age at registration')
points(h2$x, h2$y, type = "l", col = "mediumseagreen")
abline(v = median(x2, na.rm = TRUE), lty = 4, col = "mediumseagreen" )
dev.off()## quartz_off_screen
## 2x = users_demo_sympto$kg
plot_demographics_histogram(x = x, xlim = c(40,140),xlab = 'Weight (kg)')
## x y
## 1 20 4.583372e-04
## 2 21 9.166743e-05
## 3 22 2.750023e-04
## 4 23 0.000000e+00
## 5 24 9.166743e-05
## 6 25 9.166743e-05
## 7 26 2.750023e-04
## 8 27 9.166743e-05
## 9 28 9.166743e-05
## 10 29 0.000000e+00
## 11 30 9.166743e-05
## 12 31 2.750023e-04
## 13 32 2.750023e-04
## 14 33 9.166743e-05
## 15 34 2.750023e-04
## 16 35 1.833349e-04
## 17 36 0.000000e+00
## 18 37 3.666697e-04
## 19 38 4.583372e-04
## 20 39 6.416720e-04
## 21 40 1.558346e-03
## 22 41 1.375011e-03
## 23 42 2.658355e-03
## 24 43 4.766706e-03
## 25 44 4.950041e-03
## 26 45 8.800073e-03
## 27 46 1.090842e-02
## 28 47 1.448345e-02
## 29 48 2.044184e-02
## 30 49 1.723348e-02
## 31 50 4.271702e-02
## 32 51 2.218352e-02
## 33 52 4.033367e-02
## 34 53 4.225869e-02
## 35 54 3.932533e-02
## 36 55 5.115043e-02
## 37 56 3.978366e-02
## 38 57 4.528371e-02
## 39 58 4.794207e-02
## 40 59 3.492529e-02
## 41 60 6.260886e-02
## 42 61 2.245852e-02
## 43 62 3.602530e-02
## 44 63 3.474196e-02
## 45 64 2.850857e-02
## 46 65 3.602530e-02
## 47 66 1.705014e-02
## 48 67 1.677514e-02
## 49 68 2.420020e-02
## 50 69 1.173343e-02
## 51 70 2.878357e-02
## 52 71 9.625080e-03
## 53 72 1.274177e-02
## 54 73 1.182510e-02
## 55 74 7.791732e-03
## 56 75 1.860849e-02
## 57 76 6.783390e-03
## 58 77 7.150060e-03
## 59 78 8.341736e-03
## 60 79 5.041709e-03
## 61 80 1.274177e-02
## 62 81 3.391695e-03
## 63 82 5.500046e-03
## 64 83 5.133376e-03
## 65 84 3.391695e-03
## 66 85 6.875057e-03
## 67 86 4.125034e-03
## 68 87 2.566688e-03
## 69 88 3.025025e-03
## 70 89 2.016683e-03
## 71 90 6.691722e-03
## 72 91 2.016683e-03
## 73 92 1.558346e-03
## 74 93 2.200018e-03
## 75 94 1.741681e-03
## 76 95 3.758365e-03
## 77 96 2.566688e-03
## 78 97 2.108351e-03
## 79 98 1.833349e-03
## 80 99 1.283344e-03
## 81 100 3.300028e-03
## 82 101 8.250069e-04
## 83 102 1.008342e-03
## 84 103 8.250069e-04
## 85 104 4.583372e-04
## 86 105 1.283344e-03
## 87 106 5.500046e-04
## 88 107 6.416720e-04
## 89 108 9.166743e-04
## 90 109 1.191677e-03
## 91 110 1.741681e-03
## 92 111 1.833349e-04
## 93 112 2.750023e-04
## 94 113 6.416720e-04
## 95 114 7.333394e-04
## 96 115 1.375011e-03
## 97 116 3.666697e-04
## 98 117 1.833349e-04
## 99 118 4.583372e-04
## 100 119 3.666697e-04
## 101 120 1.008342e-03
## 102 121 1.833349e-04
## 103 122 1.833349e-04
## 104 123 1.833349e-04
## 105 124 2.750023e-04
## 106 125 3.666697e-04
## 107 126 0.000000e+00
## 108 127 0.000000e+00
## 109 128 1.833349e-04
## 110 129 0.000000e+00
## 111 130 2.750023e-04
## 112 131 9.166743e-05
## 113 132 0.000000e+00
## 114 133 9.166743e-05
## 115 134 2.750023e-04
## 116 135 9.166743e-05
## 117 136 1.833349e-04
## 118 137 0.000000e+00
## 119 138 1.833349e-04
## 120 139 0.000000e+00
## 121 140 9.166743e-05
## 122 141 0.000000e+00
## 123 142 0.000000e+00
## 124 143 0.000000e+00
## 125 144 0.000000e+00
## 126 145 9.166743e-05
## 127 146 0.000000e+00
## 128 147 0.000000e+00
## 129 148 0.000000e+00
## 130 149 3.666697e-04pdf.file = paste0(IO$panels,'2_demographics_weight.pdf')
pdf(pdf.file, width = 3.5, height = 2)
x = users_demo_sympto$kg
plot_demographics_histogram(x = x, xlim = c(40,140),xlab = 'Weight (kg)')## x y
## 1 20 4.583372e-04
## 2 21 9.166743e-05
## 3 22 2.750023e-04
## 4 23 0.000000e+00
## 5 24 9.166743e-05
## 6 25 9.166743e-05
## 7 26 2.750023e-04
## 8 27 9.166743e-05
## 9 28 9.166743e-05
## 10 29 0.000000e+00
## 11 30 9.166743e-05
## 12 31 2.750023e-04
## 13 32 2.750023e-04
## 14 33 9.166743e-05
## 15 34 2.750023e-04
## 16 35 1.833349e-04
## 17 36 0.000000e+00
## 18 37 3.666697e-04
## 19 38 4.583372e-04
## 20 39 6.416720e-04
## 21 40 1.558346e-03
## 22 41 1.375011e-03
## 23 42 2.658355e-03
## 24 43 4.766706e-03
## 25 44 4.950041e-03
## 26 45 8.800073e-03
## 27 46 1.090842e-02
## 28 47 1.448345e-02
## 29 48 2.044184e-02
## 30 49 1.723348e-02
## 31 50 4.271702e-02
## 32 51 2.218352e-02
## 33 52 4.033367e-02
## 34 53 4.225869e-02
## 35 54 3.932533e-02
## 36 55 5.115043e-02
## 37 56 3.978366e-02
## 38 57 4.528371e-02
## 39 58 4.794207e-02
## 40 59 3.492529e-02
## 41 60 6.260886e-02
## 42 61 2.245852e-02
## 43 62 3.602530e-02
## 44 63 3.474196e-02
## 45 64 2.850857e-02
## 46 65 3.602530e-02
## 47 66 1.705014e-02
## 48 67 1.677514e-02
## 49 68 2.420020e-02
## 50 69 1.173343e-02
## 51 70 2.878357e-02
## 52 71 9.625080e-03
## 53 72 1.274177e-02
## 54 73 1.182510e-02
## 55 74 7.791732e-03
## 56 75 1.860849e-02
## 57 76 6.783390e-03
## 58 77 7.150060e-03
## 59 78 8.341736e-03
## 60 79 5.041709e-03
## 61 80 1.274177e-02
## 62 81 3.391695e-03
## 63 82 5.500046e-03
## 64 83 5.133376e-03
## 65 84 3.391695e-03
## 66 85 6.875057e-03
## 67 86 4.125034e-03
## 68 87 2.566688e-03
## 69 88 3.025025e-03
## 70 89 2.016683e-03
## 71 90 6.691722e-03
## 72 91 2.016683e-03
## 73 92 1.558346e-03
## 74 93 2.200018e-03
## 75 94 1.741681e-03
## 76 95 3.758365e-03
## 77 96 2.566688e-03
## 78 97 2.108351e-03
## 79 98 1.833349e-03
## 80 99 1.283344e-03
## 81 100 3.300028e-03
## 82 101 8.250069e-04
## 83 102 1.008342e-03
## 84 103 8.250069e-04
## 85 104 4.583372e-04
## 86 105 1.283344e-03
## 87 106 5.500046e-04
## 88 107 6.416720e-04
## 89 108 9.166743e-04
## 90 109 1.191677e-03
## 91 110 1.741681e-03
## 92 111 1.833349e-04
## 93 112 2.750023e-04
## 94 113 6.416720e-04
## 95 114 7.333394e-04
## 96 115 1.375011e-03
## 97 116 3.666697e-04
## 98 117 1.833349e-04
## 99 118 4.583372e-04
## 100 119 3.666697e-04
## 101 120 1.008342e-03
## 102 121 1.833349e-04
## 103 122 1.833349e-04
## 104 123 1.833349e-04
## 105 124 2.750023e-04
## 106 125 3.666697e-04
## 107 126 0.000000e+00
## 108 127 0.000000e+00
## 109 128 1.833349e-04
## 110 129 0.000000e+00
## 111 130 2.750023e-04
## 112 131 9.166743e-05
## 113 132 0.000000e+00
## 114 133 9.166743e-05
## 115 134 2.750023e-04
## 116 135 9.166743e-05
## 117 136 1.833349e-04
## 118 137 0.000000e+00
## 119 138 1.833349e-04
## 120 139 0.000000e+00
## 121 140 9.166743e-05
## 122 141 0.000000e+00
## 123 142 0.000000e+00
## 124 143 0.000000e+00
## 125 144 0.000000e+00
## 126 145 9.166743e-05
## 127 146 0.000000e+00
## 128 147 0.000000e+00
## 129 148 0.000000e+00
## 130 149 3.666697e-04dev.off()## quartz_off_screen
## 2x = users_demo_sympto$cm
plot_demographics_histogram(x = x, xlim = c(140,190),xlab = 'Height (cm)')
## x y
## 1 100 9.141603e-05
## 2 101 0.000000e+00
## 3 102 0.000000e+00
## 4 103 9.141603e-05
## 5 104 0.000000e+00
## 6 105 9.141603e-05
## 7 106 0.000000e+00
## 8 107 1.828321e-04
## 9 108 0.000000e+00
## 10 109 0.000000e+00
## 11 110 0.000000e+00
## 12 111 0.000000e+00
## 13 112 0.000000e+00
## 14 113 0.000000e+00
## 15 114 0.000000e+00
## 16 115 0.000000e+00
## 17 116 0.000000e+00
## 18 117 0.000000e+00
## 19 118 0.000000e+00
## 20 119 0.000000e+00
## 21 120 2.742481e-04
## 22 121 0.000000e+00
## 23 122 0.000000e+00
## 24 123 0.000000e+00
## 25 124 0.000000e+00
## 26 125 0.000000e+00
## 27 126 0.000000e+00
## 28 127 0.000000e+00
## 29 128 0.000000e+00
## 30 129 0.000000e+00
## 31 130 9.141603e-05
## 32 131 9.141603e-05
## 33 132 0.000000e+00
## 34 133 0.000000e+00
## 35 134 0.000000e+00
## 36 135 0.000000e+00
## 37 136 0.000000e+00
## 38 137 0.000000e+00
## 39 138 0.000000e+00
## 40 139 9.141603e-05
## 41 140 2.742481e-04
## 42 141 0.000000e+00
## 43 142 0.000000e+00
## 44 143 9.141603e-05
## 45 144 3.656641e-04
## 46 145 4.570802e-04
## 47 146 8.227443e-04
## 48 147 7.313283e-04
## 49 148 1.554073e-03
## 50 149 9.141603e-04
## 51 150 8.501691e-03
## 52 151 4.205138e-03
## 53 152 6.399122e-03
## 54 153 1.197550e-02
## 55 154 9.507268e-03
## 56 155 1.608922e-02
## 57 156 1.718621e-02
## 58 157 3.062437e-02
## 59 158 3.894323e-02
## 60 159 2.788189e-02
## 61 160 7.578389e-02
## 62 161 2.742481e-02
## 63 162 5.027882e-02
## 64 163 6.399122e-02
## 65 164 5.411829e-02
## 66 165 7.715513e-02
## 67 166 3.446384e-02
## 68 167 4.954749e-02
## 69 168 7.231008e-02
## 70 169 4.954749e-02
## 71 170 7.157875e-02
## 72 171 2.541366e-02
## 73 172 3.738916e-02
## 74 173 3.245269e-02
## 75 174 2.212268e-02
## 76 175 2.248834e-02
## 77 176 1.471798e-02
## 78 177 8.501691e-03
## 79 178 1.078709e-02
## 80 179 4.479386e-03
## 81 180 8.501691e-03
## 82 181 2.193985e-03
## 83 182 2.011153e-03
## 84 183 2.285401e-03
## 85 184 8.227443e-04
## 86 185 8.227443e-04
## 87 186 4.570802e-04
## 88 187 2.742481e-04
## 89 188 1.828321e-04
## 90 189 0.000000e+00
## 91 190 9.141603e-05
## 92 191 0.000000e+00
## 93 192 0.000000e+00
## 94 193 9.141603e-05
## 95 194 0.000000e+00
## 96 195 0.000000e+00
## 97 196 9.141603e-05
## 98 197 0.000000e+00
## 99 198 0.000000e+00
## 100 199 0.000000e+00
## 101 200 0.000000e+00
## 102 201 0.000000e+00
## 103 202 0.000000e+00
## 104 203 0.000000e+00
## 105 204 0.000000e+00
## 106 205 0.000000e+00
## 107 206 0.000000e+00
## 108 207 0.000000e+00
## 109 208 0.000000e+00
## 110 209 0.000000e+00
## 111 210 0.000000e+00
## 112 211 0.000000e+00
## 113 212 0.000000e+00
## 114 213 0.000000e+00
## 115 214 0.000000e+00
## 116 215 0.000000e+00
## 117 216 0.000000e+00
## 118 217 0.000000e+00
## 119 218 0.000000e+00
## 120 219 0.000000e+00
## 121 220 0.000000e+00
## 122 221 0.000000e+00
## 123 222 0.000000e+00
## 124 223 0.000000e+00
## 125 224 0.000000e+00
## 126 225 0.000000e+00
## 127 226 0.000000e+00
## 128 227 0.000000e+00
## 129 228 0.000000e+00
## 130 229 9.141603e-05pdf.file = paste0(IO$panels,'2_demographics_height.pdf')
pdf(pdf.file, width = 3.5, height = 2)
plot_demographics_histogram(x = x, xlim = c(140,190),xlab = 'Height (cm)')## x y
## 1 100 9.141603e-05
## 2 101 0.000000e+00
## 3 102 0.000000e+00
## 4 103 9.141603e-05
## 5 104 0.000000e+00
## 6 105 9.141603e-05
## 7 106 0.000000e+00
## 8 107 1.828321e-04
## 9 108 0.000000e+00
## 10 109 0.000000e+00
## 11 110 0.000000e+00
## 12 111 0.000000e+00
## 13 112 0.000000e+00
## 14 113 0.000000e+00
## 15 114 0.000000e+00
## 16 115 0.000000e+00
## 17 116 0.000000e+00
## 18 117 0.000000e+00
## 19 118 0.000000e+00
## 20 119 0.000000e+00
## 21 120 2.742481e-04
## 22 121 0.000000e+00
## 23 122 0.000000e+00
## 24 123 0.000000e+00
## 25 124 0.000000e+00
## 26 125 0.000000e+00
## 27 126 0.000000e+00
## 28 127 0.000000e+00
## 29 128 0.000000e+00
## 30 129 0.000000e+00
## 31 130 9.141603e-05
## 32 131 9.141603e-05
## 33 132 0.000000e+00
## 34 133 0.000000e+00
## 35 134 0.000000e+00
## 36 135 0.000000e+00
## 37 136 0.000000e+00
## 38 137 0.000000e+00
## 39 138 0.000000e+00
## 40 139 9.141603e-05
## 41 140 2.742481e-04
## 42 141 0.000000e+00
## 43 142 0.000000e+00
## 44 143 9.141603e-05
## 45 144 3.656641e-04
## 46 145 4.570802e-04
## 47 146 8.227443e-04
## 48 147 7.313283e-04
## 49 148 1.554073e-03
## 50 149 9.141603e-04
## 51 150 8.501691e-03
## 52 151 4.205138e-03
## 53 152 6.399122e-03
## 54 153 1.197550e-02
## 55 154 9.507268e-03
## 56 155 1.608922e-02
## 57 156 1.718621e-02
## 58 157 3.062437e-02
## 59 158 3.894323e-02
## 60 159 2.788189e-02
## 61 160 7.578389e-02
## 62 161 2.742481e-02
## 63 162 5.027882e-02
## 64 163 6.399122e-02
## 65 164 5.411829e-02
## 66 165 7.715513e-02
## 67 166 3.446384e-02
## 68 167 4.954749e-02
## 69 168 7.231008e-02
## 70 169 4.954749e-02
## 71 170 7.157875e-02
## 72 171 2.541366e-02
## 73 172 3.738916e-02
## 74 173 3.245269e-02
## 75 174 2.212268e-02
## 76 175 2.248834e-02
## 77 176 1.471798e-02
## 78 177 8.501691e-03
## 79 178 1.078709e-02
## 80 179 4.479386e-03
## 81 180 8.501691e-03
## 82 181 2.193985e-03
## 83 182 2.011153e-03
## 84 183 2.285401e-03
## 85 184 8.227443e-04
## 86 185 8.227443e-04
## 87 186 4.570802e-04
## 88 187 2.742481e-04
## 89 188 1.828321e-04
## 90 189 0.000000e+00
## 91 190 9.141603e-05
## 92 191 0.000000e+00
## 93 192 0.000000e+00
## 94 193 9.141603e-05
## 95 194 0.000000e+00
## 96 195 0.000000e+00
## 97 196 9.141603e-05
## 98 197 0.000000e+00
## 99 198 0.000000e+00
## 100 199 0.000000e+00
## 101 200 0.000000e+00
## 102 201 0.000000e+00
## 103 202 0.000000e+00
## 104 203 0.000000e+00
## 105 204 0.000000e+00
## 106 205 0.000000e+00
## 107 206 0.000000e+00
## 108 207 0.000000e+00
## 109 208 0.000000e+00
## 110 209 0.000000e+00
## 111 210 0.000000e+00
## 112 211 0.000000e+00
## 113 212 0.000000e+00
## 114 213 0.000000e+00
## 115 214 0.000000e+00
## 116 215 0.000000e+00
## 117 216 0.000000e+00
## 118 217 0.000000e+00
## 119 218 0.000000e+00
## 120 219 0.000000e+00
## 121 220 0.000000e+00
## 122 221 0.000000e+00
## 123 222 0.000000e+00
## 124 223 0.000000e+00
## 125 224 0.000000e+00
## 126 225 0.000000e+00
## 127 226 0.000000e+00
## 128 227 0.000000e+00
## 129 228 0.000000e+00
## 130 229 9.141603e-05dev.off()## quartz_off_screen
## 2x = users_demo_sympto$bmi
plot_demographics_histogram(x = x, xlim = c(0,50),xlab = 'BMI')
## x y
## 1 8 4.592212e-04
## 2 9 4.592212e-04
## 3 10 4.592212e-04
## 4 11 3.673769e-04
## 5 12 4.592212e-04
## 6 13 3.673769e-04
## 7 14 3.673769e-04
## 8 15 1.193975e-03
## 9 16 5.969875e-03
## 10 17 2.296106e-02
## 11 18 6.300514e-02
## 12 19 1.118663e-01
## 13 20 1.379500e-01
## 14 21 1.373990e-01
## 15 22 1.221528e-01
## 16 23 8.403747e-02
## 17 24 6.998530e-02
## 18 25 4.821822e-02
## 19 26 3.655400e-02
## 20 27 3.003306e-02
## 21 28 2.406319e-02
## 22 29 1.735856e-02
## 23 30 1.726672e-02
## 24 31 9.551800e-03
## 25 32 1.092946e-02
## 26 33 8.082292e-03
## 27 34 6.520940e-03
## 28 35 6.245408e-03
## 29 36 5.510654e-03
## 30 37 5.051433e-03
## 31 38 2.663483e-03
## 32 39 2.112417e-03
## 33 40 1.377663e-03
## 34 41 2.112417e-03
## 35 42 1.561352e-03
## 36 43 1.102131e-03
## 37 44 2.755327e-04
## 38 45 9.184423e-04
## 39 46 7.347539e-04
## 40 47 5.510654e-04
## 41 48 3.673769e-04
## 42 49 9.184423e-05
## 43 50 9.184423e-05
## 44 51 9.184423e-05
## 45 52 0.000000e+00
## 46 53 2.755327e-04
## 47 54 9.184423e-05
## 48 55 9.184423e-05
## 49 56 1.836885e-04
## 50 57 1.836885e-04
## 51 58 0.000000e+00
## 52 59 0.000000e+00
## 53 60 0.000000e+00
## 54 61 0.000000e+00
## 55 62 0.000000e+00
## 56 63 0.000000e+00
## 57 64 0.000000e+00
## 58 65 0.000000e+00
## 59 66 0.000000e+00
## 60 67 9.184423e-05
## 61 68 0.000000e+00
## 62 69 9.184423e-05
## 63 70 0.000000e+00
## 64 71 0.000000e+00
## 65 72 0.000000e+00
## 66 73 0.000000e+00
## 67 74 0.000000e+00
## 68 75 0.000000e+00
## 69 76 0.000000e+00
## 70 77 0.000000e+00
## 71 78 0.000000e+00
## 72 79 0.000000e+00
## 73 80 0.000000e+00
## 74 81 0.000000e+00
## 75 82 0.000000e+00
## 76 83 0.000000e+00
## 77 84 0.000000e+00
## 78 85 0.000000e+00
## 79 86 0.000000e+00
## 80 87 0.000000e+00
## 81 88 0.000000e+00
## 82 89 0.000000e+00
## 83 90 0.000000e+00
## 84 91 0.000000e+00
## 85 92 0.000000e+00
## 86 93 0.000000e+00
## 87 94 0.000000e+00
## 88 95 0.000000e+00
## 89 96 0.000000e+00
## 90 97 0.000000e+00
## 91 98 0.000000e+00
## 92 99 0.000000e+00
## 93 100 0.000000e+00
## 94 101 0.000000e+00
## 95 102 0.000000e+00
## 96 103 0.000000e+00
## 97 104 9.184423e-05pdf.file = paste0(IO$panels,'S1C_demographics_bmi.pdf')
pdf(pdf.file, width = 3.5, height = 2)
plot_demographics_histogram(x = x, xlim = c(0,50),xlab = 'BMI')## x y
## 1 8 4.592212e-04
## 2 9 4.592212e-04
## 3 10 4.592212e-04
## 4 11 3.673769e-04
## 5 12 4.592212e-04
## 6 13 3.673769e-04
## 7 14 3.673769e-04
## 8 15 1.193975e-03
## 9 16 5.969875e-03
## 10 17 2.296106e-02
## 11 18 6.300514e-02
## 12 19 1.118663e-01
## 13 20 1.379500e-01
## 14 21 1.373990e-01
## 15 22 1.221528e-01
## 16 23 8.403747e-02
## 17 24 6.998530e-02
## 18 25 4.821822e-02
## 19 26 3.655400e-02
## 20 27 3.003306e-02
## 21 28 2.406319e-02
## 22 29 1.735856e-02
## 23 30 1.726672e-02
## 24 31 9.551800e-03
## 25 32 1.092946e-02
## 26 33 8.082292e-03
## 27 34 6.520940e-03
## 28 35 6.245408e-03
## 29 36 5.510654e-03
## 30 37 5.051433e-03
## 31 38 2.663483e-03
## 32 39 2.112417e-03
## 33 40 1.377663e-03
## 34 41 2.112417e-03
## 35 42 1.561352e-03
## 36 43 1.102131e-03
## 37 44 2.755327e-04
## 38 45 9.184423e-04
## 39 46 7.347539e-04
## 40 47 5.510654e-04
## 41 48 3.673769e-04
## 42 49 9.184423e-05
## 43 50 9.184423e-05
## 44 51 9.184423e-05
## 45 52 0.000000e+00
## 46 53 2.755327e-04
## 47 54 9.184423e-05
## 48 55 9.184423e-05
## 49 56 1.836885e-04
## 50 57 1.836885e-04
## 51 58 0.000000e+00
## 52 59 0.000000e+00
## 53 60 0.000000e+00
## 54 61 0.000000e+00
## 55 62 0.000000e+00
## 56 63 0.000000e+00
## 57 64 0.000000e+00
## 58 65 0.000000e+00
## 59 66 0.000000e+00
## 60 67 9.184423e-05
## 61 68 0.000000e+00
## 62 69 9.184423e-05
## 63 70 0.000000e+00
## 64 71 0.000000e+00
## 65 72 0.000000e+00
## 66 73 0.000000e+00
## 67 74 0.000000e+00
## 68 75 0.000000e+00
## 69 76 0.000000e+00
## 70 77 0.000000e+00
## 71 78 0.000000e+00
## 72 79 0.000000e+00
## 73 80 0.000000e+00
## 74 81 0.000000e+00
## 75 82 0.000000e+00
## 76 83 0.000000e+00
## 77 84 0.000000e+00
## 78 85 0.000000e+00
## 79 86 0.000000e+00
## 80 87 0.000000e+00
## 81 88 0.000000e+00
## 82 89 0.000000e+00
## 83 90 0.000000e+00
## 84 91 0.000000e+00
## 85 92 0.000000e+00
## 86 93 0.000000e+00
## 87 94 0.000000e+00
## 88 95 0.000000e+00
## 89 96 0.000000e+00
## 90 97 0.000000e+00
## 91 98 0.000000e+00
## 92 99 0.000000e+00
## 93 100 0.000000e+00
## 94 101 0.000000e+00
## 95 102 0.000000e+00
## 96 103 0.000000e+00
## 97 104 9.184423e-05dev.off()## quartz_off_screen
## 2rm(users_demo_sympto, x, pdf.file)3 Tracking frequency
3.1 Kindara
load(paste0(IO$kindara_data_cycles,"03_standard_CYCLES.Rdata"), verbose = TRUE)## Loading objects:
## CYCLESCYCLES = CYCLES[(CYCLES$tracking_freq>=0) & (CYCLES$tracking_freq <= 1 ), ]get_hist <- function(p) {
d <- ggplot_build(p)$data[[1]]
data.frame(x = d$x, xmin = d$xmin, xmax = d$xmax, y = d$y, group = d$group)
}sex_activity = c("no logged sex","logged sex")
CYCLES$has_logged_sex = factor(sex_activity[(CYCLES$n_sex_tot > 0)+1], levels = sex_activity)
a = aggregate(tracking_freq ~ has_logged_sex, CYCLES, function(x) round(median(x),digits = 2))
g = ggplot(CYCLES, aes(x = tracking_freq, fill = has_logged_sex))+ #coord_flip()+
geom_histogram(binwidth = 0.025,aes(y = stat(width*density))) +
geom_vline(data = a, aes(xintercept = tracking_freq, group = has_logged_sex), linetype = 2)+
geom_text(data = a, aes(x = tracking_freq, y = 0.25, label = tracking_freq , group = has_logged_sex), hjust = -0.5)+
facet_grid(. ~ has_logged_sex) +
expand_limits(x = 0)+
scale_x_continuous(breaks = seq(0,1,by = 0.2))+
scale_y_continuous(labels = percent_format()) +
xlab("Tracking frequency (#day with obs. / cycle length)")+
ylab("% of cycles")+
guides(fill = FALSE)+
theme_light()
g
ggsave(filename = paste0(IO$panels,"2_tracking_frequency_kindara.pdf"),width = 4.5, height = 3)
h_values = get_hist(g)
h_values$group_name = sex_activity[h_values$group]
h_values$tracking_freq = h_values$x
h_values$perc_of_cycles = round(h_values$y * 100, digits = 2)
save(h_values, file = paste0(IO$output_Rdata, "tracking_frequency_kindara.Rdata"))
write.csv(h_values[,c("group_name","tracking_freq","perc_of_cycles")], file = paste0(IO$output_csv, "tracking_frequency_kindara.csv"), row.names = FALSE)
rm(g, sex_activity, h_values, a, CYCLES)3.2 Sympto
source("Scripts/_setup_Sympto.R")load(paste0(IO$sympto_data_02_standard_cycles,'cycles.Rdata'), verbose = TRUE)## Loading objects:
## cyclesbreaks = seq(0,1,by = 0.025)
N = length(breaks)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
for(goal in 1:4){
for(sex in c(TRUE, FALSE)){
tfh = plot_tracking_freq(cycles, goal, sex)
save(tfh, file = paste0(IO$output_Rdata,'tracking_frequency_sympto_',goal,'_',sex,'.Rdata'))
write.csv(tfh, file = paste0(IO$output_csv,'tracking_frequency_sympto_',goal,'_',sex,'.csv'))
pdf.file = paste0(IO$panels,'2_tracking_frequency_sympto_',goal,'_',sex,'.pdf')
pdf(pdf.file, width = 4.5, height = 1.5)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
plot_tracking_freq(cycles, goal, sex)
dev.off()
}
}
4 Observation profiles
4.1 Examples
4.1.1 Kindara
source("Scripts/_setup_Kindara.R")load(paste0(IO$kindara_data_days,"03 selected/day10-11.Rdata"), verbose = TRUE)## Loading objects:
## dayscycle_ids = c('c3eeacfe8faf436db8bbc5972c25d208_18','c3eeacfe8faf436db8bbc5972c25d208_19','b4bf193f9ab94acead43a0eb1f6a0dcc_3','d71e51e3114b4b5fb9f7124115c901dd_4')
for(cycle_id in cycle_ids){
k = which(days$cycle_id == cycle_id)
plot_cycle(subtable = days[k,], show_temp_time = FALSE, show_weights = TRUE, show_goal = FALSE)
}
cycle_ids = c('c3eeacfe8faf436db8bbc5972c25d208_19')
for(cycle_id in cycle_ids){
k = which(days$cycle_id == cycle_id)
pdf(paste0(IO$panels,'3_cycle_example',cycle_id,'.pdf'),
useDingbats = FALSE, width = viz$pdf.full.width/1.5, height = viz$pdf.full.width/3)
plot_cycle(subtable = days[k,], show_temp_time = FALSE, show_weights = TRUE, show_goal = FALSE)
dev.off()
}
rm(cycle_ids, cycle_id, k, days)4.1.2 Sympto
source("Scripts/_setup_Sympto.R")load(paste0(IO$sympto_data_02_standard_cycles,"cycledays.Rdata"), verbose = TRUE)## Loading objects:
## cycledaysuser_id_str = c("S01370x1jgg","S00278x2rao")
cycle_nb = c(5,66)
cycle_ids = paste0(cycledays$user_id[match( user_id_str,cycledays$user_id_str)], "_",cycle_nb)
for(cycle_id in cycle_ids){
k = which(cycledays$cycle_id == cycle_id)
plot_cycle(cycledays[k,])
pdf(paste0(IO$panels,'3_cycle_example',cycle_id,'.pdf'),
useDingbats = FALSE, width = viz$pdf.full.width/1.5, height = viz$pdf.full.width/3)
plot_cycle(cycledays[k,])
dev.off()
}
rm(user_id_str,cycle_nb,cycle_ids, cycle_id, k, cycledays)4.2 Temperature
4.2.1 Kindara
source("Scripts/_setup_Kindara.R")load(file = paste0(IO$output_Rdata,"delta_BBT_distribution_kindara.Rdata"), verbose = TRUE)## Loading objects:
## delta_BBT_distribution_kindaraload(file = paste0(IO$output_Rdata,"delta_BBT_summary_kindara.Rdata"), verbose = TRUE)## Loading objects:
## delta_BBT_summary_kindaramids = as.numeric(colnames(delta_BBT_distribution_kindara[-1]))
cols = paste0('gray',as.numeric(100-round(100*as.matrix(delta_BBT_distribution_kindara[,-1])/max(delta_BBT_distribution_kindara[,-1]))))
cycleday_from_end = delta_BBT_distribution_kindara$cycleday
pdf(file = paste0(IO$panels, "3_temperature_profile_points_end_kindara.pdf"), width = 7, height = 3, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
plot(rep.col(delta_BBT_distribution_kindara$cycleday, ncol(delta_BBT_distribution_kindara)-1),
rep.row(mids,nrow(delta_BBT_distribution_kindara)),
pch = 16,
col = cols,
xlab= 'Days from next cycle',
ylab= expression(paste(Delta,'', T,' from T-pc25 [F] ')),
axes = FALSE)
axis(1, at = c(-1,seq(-7,-31,by= -7)))
axis(2)
axis(4, at = round(range(delta_BBT_summary_kindara$median),digits = 2))
box()
#10-90
polygon(c(cycleday_from_end,rev(cycleday_from_end)),
c(delta_BBT_summary_kindara$p.10, rev(delta_BBT_summary_kindara$p.90)),border = NA, col = rgb(0,0.5,0.8,0.25))
#25-75
polygon(c(cycleday_from_end,rev(cycleday_from_end)),
c(delta_BBT_summary_kindara$p.25, rev(delta_BBT_summary_kindara$p.75)),border = NA, col = rgb(0,0.5,0.8,0.35))
# median
points(cycleday_from_end, delta_BBT_summary_kindara$median, type = 'l', lwd = 3, col = rgb(0,0.5,0.8))
segments(x0 = cycleday_from_end[which.min(delta_BBT_summary_kindara$median)],
y0 = min(delta_BBT_summary_kindara$median),
x1 = 0, y1 = min(delta_BBT_summary_kindara$median), col = 'gray10', lty = 2)
segments(x0 = cycleday_from_end[which.max(delta_BBT_summary_kindara$median)],
y0 = max(delta_BBT_summary_kindara$median),
x1 = 0, y1 = max(delta_BBT_summary_kindara$median), col = 'gray10', lty = 2)
dev.off()## quartz_off_screen
## 2rm(delta_BBT_distribution_kindara, delta_BBT_summary_kindara, mids, cols, cycleday_from_end)4.2.2 Sympto
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$output_Rdata,"delta_BBT_distribution_sympto.Rdata"), verbose = TRUE)## Loading objects:
## delta_BBT_distribution_symptoload(file = paste0(IO$output_Rdata,"delta_BBT_summary_sympto.Rdata"), verbose = TRUE)## Loading objects:
## delta_BBT_summary_symptopdf(file = paste0(IO$panels, "3_temperature_profile_points_end_Sympto.pdf"), width = 7, height = 3, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
plot(delta_BBT_distribution_sympto$cycleday_from_end, delta_BBT_distribution_sympto$low_norm_temp,
#xlim = c(-32,-1),
ylim = c(-0.2,1),
xlab= 'Days from next cycle',
ylab= expression(paste(Delta,'', T^o,' from T - pc25 ')),
axes = FALSE, type = 'n')
points(delta_BBT_distribution_sympto$cycleday_from_end,delta_BBT_distribution_sympto$low_norm_temp,
col = rgb(0,0,0,0.7/max(delta_BBT_distribution_sympto$density)*delta_BBT_distribution_sympto$density),cex = 1, pch = 16)
axis(1, at = c(-1,seq(-7,-31,by= -7)))
axis(2)
axis(4, at = round(range(delta_BBT_summary_sympto$median),digits = 2))
box()
#10-90
polygon(c(delta_BBT_summary_sympto$cycleday, rev(delta_BBT_summary_sympto$cycleday)),
c(delta_BBT_summary_sympto$p.10, rev(delta_BBT_summary_sympto$p.90)),border = NA, col = rgb(0,0.5,0.8,0.25))
#25-75
polygon(c(delta_BBT_summary_sympto$cycleday,rev(delta_BBT_summary_sympto$cycleday)),
c(delta_BBT_summary_sympto$p.25, rev(delta_BBT_summary_sympto$p.75)),border = NA, col = rgb(0,0.5,0.8,0.35))
# median
points(delta_BBT_summary_sympto$cycleday, delta_BBT_summary_sympto$median, type = 'l', lwd = 3, col = rgb(0,0.5,0.8))
segments(x0 = delta_BBT_summary_sympto$cycleday[which.min(delta_BBT_summary_sympto$median)],
y0 = min(delta_BBT_summary_sympto$median),
x1 = 0, y1 = min(delta_BBT_summary_sympto$median), col = 'gray10', lty = 2)
segments(x0 = delta_BBT_summary_sympto$cycleday[which.max(delta_BBT_summary_sympto$median)],
y0 = max(delta_BBT_summary_sympto$median),
x1 = 0, y1 = max(delta_BBT_summary_sympto$median), col = 'gray10', lty = 2)
dev.off()## quartz_off_screen
## 2rm(delta_BBT_distribution_sympto, delta_BBT_summary_sympto)4.2.3 Temperature by orifice (Sympto only)
load(paste0(IO$sympto_data_02_standard_cycles_with_HMM_res,"cycledays.Rdata"), verbose = TRUE)## Loading objects:
## cycledaysload(paste0(IO$sympto_data_02_standard_cycles_with_HMM_res,"cycles.Rdata"), verbose = TRUE)## Loading objects:
## cycles######
# temperature distribution according to orifice
m = match(cycledays$cycle_id, cycles$cycle_id)
cycledays$orifice = cycles$temp_method_txt[m]
orifice = factor(cycledays$orifice, unique(cycledays$orifice))
phase = factor(cycledays$out_phase, levels = c(1,2,4,3))
pdf(file = paste0(IO$panels, "S2_temperature_by_orifice.pdf"), width = 6, height = 3.5, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
b = boxplot(cycledays$temp_c ~ orifice, ylim = c(35.8, 38),
col = viz_cols$colors_lsy[c(2,4,8,10)], pch = 16, border = viz_cols$colors_lsy_dark[c(2,4,8,10)], cex = 0.5)
dev.off()## quartz_off_screen
## 2pdf(file = paste0(IO$panels, "S2_temperature_by_orifice_and_phase.pdf"), width = 12, height = 3.5, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
b = boxplot(cycledays$temp_c[orifice != "Undefined"] ~ phase[orifice != "Undefined"]:orifice[orifice != "Undefined"], ylim = c(35.8, 38),
col = dict$phases$colors[match(levels(phase),dict$phases$index) ], pch = 16, border = 'black', cex = 0.5)
dev.off()## quartz_off_screen
## 2rm(m, orifice, phase, cycledays, cycles)4.3 Bleeding
4.3.1 Kindara
source("Scripts/_setup_Kindara.R")load(file = paste0(IO$output_Rdata,"bleeding_distribution_kindara.Rdata"), verbose = TRUE)## Loading objects:
## bleeding_distribution_kindaraat.y = seq(0,1,by = 0.25)
pdf(file = paste0(IO$panels, "3_bleeding_profile_kindara.pdf"), width = 8, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 3) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
b = barplot(t(as.matrix(bleeding_distribution_kindara[,-1][,c(2:5,1)])),
col = dict$bleeding$colors[c(2:5,1)], border = NA,
xlab = 'Days in cycle', ylab = '',
axes = FALSE,ylim = c(0,1),
names.arg = bleeding_distribution_kindara$cycleday, add = FALSE, plot = TRUE)
axis(2, at = seq(0,1,by = 0.25), labels = seq(0,100,by = 25))
axis(4, at = seq(0,1,by = 0.25), labels = seq(0,100,by = 25))
dev.off()## quartz_off_screen
## 2rm(bleeding_distribution_kindara, at.y, b)4.3.2 Sympto
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$output_Rdata,"bleeding_distribution_sympto.Rdata"), verbose = TRUE)## Loading objects:
## bleeding_distribution_symptopdf(file = paste0(IO$panels, "3_bleeding_profile_Sympto.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
b = barplot(t(bleeding_distribution_sympto[,-1]),
col = dict$bleeding$colors[2:4], border = NA,
xlab = 'Days in cycle', ylab = 'Frequency',
axes = FALSE,ylim = c(0,1),
names.arg = bleeding_distribution_sympto$cycleday, add = FALSE, plot = TRUE)
axis(2, at = seq(0,1,by = 0.25))
dev.off()## quartz_off_screen
## 2rm(bleeding_distribution_sympto)4.4 Mucus
4.4.1 Kindara
source("Scripts/_setup_Kindara.R")load(file = paste0(IO$output_Rdata,"cervical_mucus_distribution_kindara.Rdata"), verbose = TRUE)## Loading objects:
## cervical_mucus_distribution_kindaraat.y = seq(0,1, by = 0.025)
pdf(file = paste0(IO$panels, "3_mucus_profile_end_lines_kindara.pdf"), width = 7, height = 2.5, useDingbats = FALSE)
par(las = 1, mar = c(3, 2.5, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
matplot(cervical_mucus_distribution_kindara[,-1], type = 'l',
col = dict$mucus$colors[which(dict$mucus$index>=0)],
lwd = 1.5*dict$mucus$cex[which(dict$mucus$index>=0)],
lty = dict$mucus$lty[which(dict$mucus$index>=0)],
xlab = 'Days from next cycle', ylab = '',
axes = FALSE,ylim = c(0,0.10))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(1,22,by = 7),28), labels = c(seq(-28,-7,by=7),1))
dev.off()## quartz_off_screen
## 2rm(cervical_mucus_distribution_kindara, at.y)4.4.2 Sympto
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$output_Rdata,"cervical_mucus_distribution_sympto.Rdata"), verbose = TRUE)## Loading objects:
## cervical_mucus_distribution_symptopdf(file = paste0(IO$panels, "3_mucus_profile_Sympto.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
at.y = seq(0,1,by = 0.25)
matplot(x = cervical_mucus_distribution_sympto$cycleday,
y = cervical_mucus_distribution_sympto[,-1],
type = 'l', lwd = 2, lty = 1,
col = dict$mucus$colors[match(colnames(cervical_mucus_distribution_sympto[-1]),dict$mucus$hmm.symbols)],
xlab = 'Days in cycle', ylab = '',
axes = FALSE,ylim = c(0,0.5))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(-28,-7,by = 7),-1))
dev.off()## quartz_off_screen
## 2rm(at.y,cervical_mucus_distribution_sympto)4.5 Cervix
4.5.1 Kindara
source("Scripts/_setup_Kindara.R")load(file = paste0(IO$output_Rdata,"cervix_distribution_kindara.Rdata"), verbose = TRUE)## Loading objects:
## cervix_distribution_kindaraat.y = seq(0,1,by = 0.005)
# openness
pdf(file = paste0(IO$panels, "S2_cervix_openness_profile_end_lines.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 2.5, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
matplot(cervix_distribution_kindara[,grep("o_",colnames(cervix_distribution_kindara))], type = 'l', lwd = 2, lty = 1,
col = c(viz_cols$yellow, viz_cols$light.blue, viz_cols$blue),
xlab = 'Days from next cycle', ylab = '',
axes = FALSE,ylim = c(0,0.025))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(1,22,by = 7),28), labels = c(seq(-28,-7,by=7),1))
dev.off()## quartz_off_screen
## 2# height
pdf(file = paste0(IO$panels, "S2_cervix_height_profile_end_lines.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 2.5, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
matplot(cervix_distribution_kindara[,grep("h_",colnames(cervix_distribution_kindara))], type = 'l', lwd = 2, lty = 1,
col = c(viz_cols$yellow, viz_cols$light.blue, viz_cols$blue),
xlab = 'Days from next cycle', ylab = '',
axes = FALSE,ylim = c(0,0.025))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(1,22,by = 7),28), labels = c(seq(-28,-7,by=7),1))
dev.off()## quartz_off_screen
## 2# firmness
pdf(file = paste0(IO$panels, "S2_cervix_firmness_profile_end_lines.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 2.5, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
matplot(cervix_distribution_kindara[,grep("f_",colnames(cervix_distribution_kindara))], type = 'l', lwd = 2, lty = 1,
col = c(viz_cols$yellow, viz_cols$light.blue, viz_cols$blue),
xlab = 'Days from next cycle', ylab = '',
axes = FALSE,ylim = c(0,0.025))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(1,22,by = 7),28), labels = c(seq(-28,-7,by=7),1))
dev.off()## quartz_off_screen
## 2rm(cervix_distribution_kindara, at.y)
4.5.2 Sympto
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$output_Rdata,"cervix_distribution_sympto.Rdata"), verbose = TRUE)## Loading objects:
## cervix_distribution_symptopdf(file = paste0(IO$panels, "S2_cervix_profile_Sympto.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
at.y = seq(0, 1, by = 0.01)
matplot(x = cervix_distribution_sympto$cycleday,
y = cervix_distribution_sympto[,-1],
type = 'l', lwd = 2, lty = 1,
col = dict$cervix$colors[match(colnames(cervix_distribution_sympto[-1]),dict$cervix$names)],
xlab = 'Days in cycle', ylab = '',
axes = FALSE,ylim = c(0,0.05))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(-28,-7,by = 7),-1))
dev.off()## quartz_off_screen
## 2rm(at.y,cervix_distribution_sympto)4.6 Vaginal sensation
4.6.1 Kindara
source("Scripts/_setup_Kindara.R")load(file = paste0(IO$output_Rdata,"vaginal_sensation_distribution_kindara.Rdata"), verbose = TRUE)## Loading objects:
## vaginal_sensation_distribution_kindaraat.y = seq(0,1,by = 0.00005)
pdf(file = paste0(IO$panels, "S2_vaginal_sensation_profile_end_lines_kindara.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3.5, 0.25, 3.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
matplot(vaginal_sensation_distribution_kindara[,-1], type = 'l', lwd = 2, lty = 1,
col = c(viz_cols$light.creamy,viz_cols$yellow, viz_cols$light.blue, viz_cols$blue),
xlab = 'Days from next cycle', ylab = '',
axes = FALSE,ylim = c(0,0.0002))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(1,22,by = 7),28), labels = c(seq(-28,-7,by=7),1))
dev.off()## quartz_off_screen
## 2rm(vaginal_sensation_distribution_kindara, at.y)4.6.2 Sympto
source("Scripts/_setup_Sympto.R")load(file = paste0(IO$output_Rdata,"vaginal_sensation_distribution_sympto.Rdata"), verbose = TRUE)## Loading objects:
## vaginal_sensation_distribution_symptopdf(file = paste0(IO$panels, "S2_vaginal_sensation_profile_Sympto.pdf"), width = 7, height = 2, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.8,0.6,0), tcl = -0.4)
at.y = seq(0, 1, by = 0.05)
matplot(x = vaginal_sensation_distribution_sympto$cycleday,
y = vaginal_sensation_distribution_sympto[,-1],
type = 'l', lwd = 2, lty = 1,
col = dict$cervix$colors[match(colnames(vaginal_sensation_distribution_sympto[-1]),dict$cervix$names)],
xlab = 'Days in cycle', ylab = '',
axes = FALSE,ylim = c(0,0.25))
axis(2, at = at.y, labels = paste0(100*at.y, '%'))
axis(4, at = at.y, labels = paste0(100*at.y, '%'))
axis(1, at = c(seq(-28,-7,by = 7),-1))
dev.off()## quartz_off_screen
## 2rm(at.y,vaginal_sensation_distribution_sympto)5 Ovulation Estimation Results
5.1 Examples
5.1.1 Kindara
source("Scripts/_setup_Kindara.R")load(paste0(IO$kindara_data_days,'/04 HMM/day10-11.Rdata'), verbose = TRUE)## Loading objects:
## daysload(paste0(IO$kindara_data_cycles,'05 standard with ovu est CYCLES.Rdata'), verbose = TRUE)## Loading objects:
## CYCLEScycle_ids = c('1497e11dfe844e86b2c84ec11f557af3_4',
'201790bb7b6b42d89d05e37f656087c4_14',
'4daf807bfb9e44fc8f75e14291a844a3_30',
'8067a42a1f894de7b6f04d5fab2b6cc4_9',
'9a0f0f797020487aa14c75b9d76e6ac1_5',
'a7a672481c594ed48897caa3ce9fdb46_3',
'a7be6f1d198d406e8542a3fa6f3e7710_7',
'bb42baba52d84b458520aaf634c908fa_2',
'd4dfe0a4ee2948c59ae19f2d0b985a0f_5',
'd71e51e3114b4b5fb9f7124115c901dd_4',
'e54d03622e444b3bb5c15813bb9a2bf9_20',
'eb8ac1b2ffbd448f92f410560df616fb_3',
'fc38e35706474f87816d4f31c1ff78c8_29'
)
for(cycle_id in cycle_ids){
k = which(days$cycle_id == cycle_id)
plot_fit_results_hmm(cycletable = days[k,], cycles = CYCLES, hmm.res = NA, show_temp_time = FALSE)
pdf(paste0(IO$panels,'4_kindara_fit_cycle_example_',cycle_id,'.pdf'),
useDingbats = FALSE, width = viz$pdf.full.width/2, height = 5)
plot_fit_results_hmm(cycletable = days[k,], cycles = CYCLES, hmm.res = NA, show_temp_time = FALSE)
dev.off()
}## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots## before matplots
## before matplots5.1.2 Sympto
source("Scripts/_setup_Sympto.R")load(paste0(IO$sympto_data_02_standard_cycles_with_HMM_res,'cycledays.Rdata'), verbose = TRUE)## Loading objects:
## cycledaysload(paste0(IO$sympto_data_02_standard_cycles_with_HMM_res,'cycles.Rdata'), verbose = TRUE)## Loading objects:
## cyclescycle_ids = c('514_3','5017_22',"7197_3","990_17","7197_2")
for(cycle_id in cycle_ids){
cat(cycle_id,"\n")
j = which(cycledays$cycle_id == cycle_id)
cycletable = cycledays[j, ]
N = nrow(cycletable)
cat(N,"\n")
res = most_likely_day_of_ovulation_hmm(cycletable = cycletable)
plot_fit_results_hmm(cycletable = cycletable, hmm.res = res)
pdf(paste0(IO$panels,'4_sympto_fit_cycle_example_',cycle_id,'.pdf'),
useDingbats = FALSE, width = viz$pdf.full.width/2, height = 5)
plot_fit_results_hmm(cycletable = cycletable, hmm.res = res)
dev.off()
}## 514_3
## 28
## hM hM lM lM lM lE lE lE lE lE hE hE hE O Rise hP hP hP hP hP hP hP hP hP hP hP hP lP end## before matplots
## before matplots
## 5017_22
## 34
## hM hM hM lM lM lM lM lM lM lM lE lE lE lE lE hE hE hE hE hE O hP hP hP lP lP lP lP lP lP lP lP lP lP end## before matplots
## before matplots
## 7197_3
## 28
## hM hM hM lM lM lE lE lE lE lE lE lE lE hE O hP hP hP hP hP hP lP lP lP lP lP lP lP end## before matplots
## before matplots
## 990_17
## 25
## hM lM lE lE lE lE lE hE lE hE lE hE O hP hP hP hP hP hP hP hP hP lP lP lP end## before matplots
## before matplots
## 7197_2
## 38
## hM hM hM lM lM lE lE lE lE lE lE lE lE lE lE lE lE lE lE lE lE lE hE O hP hP hP hP hP hP hP hP lP lP lP lP lP lP end## before matplots
## before matplotsrm(cycle_id, cycle_ids, j, N, cycletable,res, cycles, cycledays)5.2 Cycle phases duration
load(file = paste0(IO$output_Rdata, 'cycle_phases_duration_distribution_kindara.Rdata') , verbose = TRUE)## Loading objects:
## cycle_phases_duration_distribution_kindarahk = cycle_phases_duration_distribution_kindaraload(file = paste0(IO$output_Rdata, 'cycle_phases_duration_distribution_sympto.Rdata') , verbose = TRUE)## Loading objects:
## cycle_phases_duration_distribution_symptohs = cycle_phases_duration_distribution_sympto
load(paste0(IO$sympto_data_03_reliable_ovulation_estimation,"cycles.Rdata"), verbose = TRUE)## Loading objects:
## cyclespdf.height = 3
pdf.width = 7
at.y = seq(0,0.22,by = 0.02)
at.x = seq(0,65,by = 7)
# cycle length
pdf(file = paste0(IO$panels, "4_cycle_length_distribution.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$cycle_length, bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Cycle length',ylab = '% of cycles',
col = 'tomato2', border = 'tomato3', bars = FALSE)
points(x = hk$duration, y = hk$cycle_length, col = 'tomato3', lw = 2, type = "l" )
points(x = hk$duration, y = hk$cycle_length, col = 'tomato3', cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('Kindara','Sympto'),
col = "tomato3",
lwd = c(2,1),lty = 1,
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2pdf(file = paste0(IO$panels, "4_estimated_ovulation_time_distribution.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$ovu, bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Estimated ovulation',ylab = '% of cycles',
col = 'gray30', border = 'gray25', bars = FALSE)
points(x = hk$duration, y = hk$estimated_ovulation, col = 'gray25', lw = 2, type = "l" )
points(x = hk$duration, y = hk$estimated_ovulation, col = 'gray25', cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('Kindara','Sympto'),
col = "gray25",
lwd = c(2,1),lty = 1,
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2pdf(file = paste0(IO$panels, "4_estimated_luteal_phase_duration_distribution.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$luteal.duration, bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Estimated luteal phase duration',ylab = '% of cycles',
col = viz_cols$yellow, border = viz_cols$dark.yellow, bars = FALSE)
points(x = hk$duration, y = hk$estimated_luteal_phase, col = viz_cols$dark.yellow, lw = 2, type = "l" )
points(x = hk$duration, y = hk$estimated_luteal_phase, col = viz_cols$dark.yellow, cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('Kindara','Sympto'),
col = viz_cols$dark.yellow,
lwd = c(2,1),lty = 1,
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2rm(pdf.height, pdf.width, at.y, at.x)pdf.height = 3
pdf.width = 7
at.y = seq(0,0.22,by = 0.02)
at.x = seq(0,65,by = 7)
# cycle length
pdf(file = paste0(IO$panels, "4_cycle_length_distribution_sympto_goals.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$cycle_length[cycles$goal_txt == "Conception"], bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Cycle length',ylab = '% of cycles',
col = 'tomato2', border = 'tomato3', bars = FALSE, plot.median = FALSE)
points(x = hs$duration, y = hs$cycle_length_avoid_preg, col = 'tomato3', lty = 2, type = "l" )
points(x = hs$duration, y = hs$cycle_length_avoid_preg, col = 'tomato3', cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('contraception','conception'),
col = 'tomato3',
lwd = 1,lty = c(2,1),
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2pdf(file = paste0(IO$panels, "4_estimated_ovulation_time_distribution_sympto_goals.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$ovu[cycles$goal_txt == "Conception"], bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Estimated ovulation',ylab = '% of cycles',
col = 'gray30', border = 'gray25', bars = FALSE)
points(x = hs$duration, y = hs$estimated_ovulation_avoid_preg, col = 'gray25', lty = 2, type = "l" )
points(x = hs$duration, y = hs$estimated_ovulation_avoid_preg, col = 'gray25', cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('contraception','conception'),
col = 'gray25',
lwd = 1,lty = c(2,1),
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2pdf(file = paste0(IO$panels, "4_estimated_luteal_phase_duration_distribution_sympto_goals.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
h = plot_hist(x = cycles$luteal.duration[cycles$goal_txt == "Conception"], bin.size = 1,
at.x = seq(0,65,by = 7), at.y = at.y,
xlab = 'Estimated luteal phase duration',ylab = '% of cycles',
col = viz_cols$yellow, border = viz_cols$dark.yellow, bars = FALSE)
points(x = hs$duration, y = hs$estimated_luteal_phase_avoid_preg, col = viz_cols$dark.yellow, lty = 2, type = "l" )
points(x = hs$duration, y = hs$estimated_luteal_phase_avoid_preg, col = viz_cols$dark.yellow, cex = 0.5, pch = 16, type = "p" )
legend('topright', legend = c('contraception','conception'),
col = viz_cols$dark.yellow,
lwd = 1,lty = c(2,1),
bty = 'o', bg ='white', box.col = 'lightgray' )
dev.off()## quartz_off_screen
## 2rm(pdf.height, pdf.width, at.y, at.x)6 Estimation parameters: Temperature shift, uncertainty and confidence score
pdf.height = 2
pdf.width = 3.56.1 Kindara
load(file = paste0(IO$output_Rdata, 'temperature_shift_distribution_kindara.Rdata'), verbose = TRUE)## Loading objects:
## temperature_shift_distribution_kindaraDT = temperature_shift_distribution_kindara
load(file = paste0(IO$output_Rdata, 'uncertainties_distribution_kindara.Rdata'), verbose = TRUE)## Loading objects:
## uncertainties_distribution_kindarasd = uncertainties_distribution_kindara
load(file = paste0(IO$output_Rdata, 'confidence_scores_distribution_kindara.Rdata'), verbose = TRUE)## Loading objects:
## confidence_scores_distribution_kindaracs = confidence_scores_distribution_kindaraat.x = DT$temperature_shifts
at.y = seq(0,1, by = 0.1)
mids = DT$temperature_shifts
by = mean(diff(mids))
pdf(file = paste0(IO$panels, "S4_DT_distribution_all_kindara.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
f = DT$fraction_of_cycles
N = length(f)
plot(mids,f, type = 'n',
axes = FALSE,xlim = range(mids-by/2, mids+by/2), ylim = c(0,0.6),
yaxs = 'i', xaxs = 'i',
xlab = 'Temperature shifts [F]', ylab = '% of cycles')
axis(1, at = at.x , col = 'lightgray', col.ticks = 'lightgray', las = 1)
axis(2, at = at.y , labels = 100*at.y, col = 'lightgray', col.ticks = 'lightgray', las = 1)
abline(h = at.y, col = 'lightgray')
abline(v = at.x, col = 'lightgray')
rect(xright =mids-by/2,ybottom = rep(0,N),xleft = mids+by/2, ytop = f,col = 'gray50', border = 'gray40')
f.0 = f; f.0[1] = 0
x.cumsum = cumsum(f.0)
med = mids[which(x.cumsum > (sum(f.0)/2))[1]]
abline(v = med, lty = 2)
text(x = med, y = at.y[2]/3, pos = 4, labels = round(med,digits = 2))
dev.off()## quartz_off_screen
## 2mids = sd$uncertainties
by = mean(diff(mids))
at.x = seq(0,max(mids), by = 1)
at.y = seq(0,1, by = 0.01)
pdf(file = paste0(IO$panels, "S4_uncertainty_distribution_all_kindara.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
f = sd$fraction_of_cycles
N = length(f)
plot(mids,f, type = 'n',
axes = FALSE,xlim = range(mids-by/2, mids+by/2), ylim = c(0,0.06),
yaxs = 'i', xaxs = 'i',
xlab = 'Uncertainty on ovu. estimation (in days)', ylab = '% of cycles')
axis(1, at = at.x , col = 'lightgray', col.ticks = 'lightgray', las = 1)
axis(2, at = at.y , labels = 100*at.y, col = 'lightgray', col.ticks = 'lightgray', las = 1)
abline(h = at.y, col = 'lightgray')
abline(v = at.x, col = 'lightgray')
rect(xright =mids-by/2,ybottom = rep(0,N),xleft = mids+by/2, ytop = f,col = 'gray50', border = 'gray40')
f.0 = f; f.0[1] = 0
x.cumsum = cumsum(f.0)
med = mids[which(x.cumsum > (sum(f.0)/2))[1]]
abline(v = med, lty = 2)
text(x = med, y = at.y[2]/3, pos = 4, labels = round(med,digits = 2))
dev.off()## quartz_off_screen
## 2mids = cs$confidence_score
by = mean(diff(mids))
at.x = seq(0,max(mids), by = 0.1)
at.y = seq(0,1, by = 0.02)
pdf(file = paste0(IO$panels, "S4_confidence_score_distribution_all_kindara.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
f = cs$fraction_of_cycles
N = length(f)
plot(mids,f, type = 'n',
axes = FALSE,xlim = range(mids-by/2, mids+by/2), ylim = c(0,0.16),
yaxs = 'i', xaxs = 'i',
xlab = 'Confidence score on ovu. estimates', ylab = '% of cycles')
axis(1, at = at.x , col = 'lightgray', col.ticks = 'lightgray', las = 1)
axis(2, at = at.y , labels = 100*at.y, col = 'lightgray', col.ticks = 'lightgray', las = 1)
abline(h = at.y, col = 'lightgray')
abline(v = at.x, col = 'lightgray')
rect(xright =mids-by/2,ybottom = rep(0,N),xleft = mids+by/2, ytop = f,col = 'gray50', border = 'gray40')
f.0 = f; f.0[1] = 0
x.cumsum = cumsum(f.0)
med = mids[which(x.cumsum > (sum(f.0)/2))[1]]
abline(v = med, lty = 2)
text(x = med, y = at.y[2]/3, pos = 4, labels = round(med,digits = 2))
dev.off()## quartz_off_screen
## 26.2 Sympto
source("Scripts/_setup_Sympto.R")pdf.height = 2
pdf.width = 3.5file = paste0(IO$sympto_data_02_standard_cycles_with_HMM_res,"cycles.Rdata")
load(file, verbose = TRUE)## Loading objects:
## cycles# DELTA T
########################
pdf(file = paste0(IO$panels, "S4_DT_distribution_all_sympto.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
x = cycles$D.T[cycles$D.T >= 0]
h.DT = plot_hist(x = x, bin.size = 0.05,
at.x = seq(0,1,by = 0.1), at.y = seq(0,0.25,by = 0.05),
xlab = 'Temperature shifts',ylab = '% of cycles' )
dev.off()## quartz_off_screen
## 2# SD on OVU ESTIMATION
########################
pdf(file = paste0(IO$panels, "S4_uncertainty_distribution_all_sympto.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
x = cycles$ovu.sd[cycles$ovu.sd >= 0]
h.sd = plot_hist(x = x, bin.size = 0.05,
at.x = seq(0,4,by = 0.5), at.y = seq(0,0.1,by = 0.02),
xlab = 'Uncertainty on ovu. estimates (in days)',ylab = '% of cycles' )
dev.off()## quartz_off_screen
## 2# CONFIDENCE
########################
pdf(file = paste0(IO$panels, "S4_confidence_score_distribution_all_sympto.pdf"), width = pdf.width, height = pdf.height, useDingbats = FALSE)
par(las = 1, mar = c(3, 3, 0.25, 2.5) + 0.1, mgp = c(1.6,0.6,0), tcl = -0.4, col.axis = 'gray60', col.lab = 'gray40')
x = cycles$confidence[cycles$confidence >= 0]
h.c = plot_hist(x = x, bin.size = 0.05,
at.x = seq(0,1,by = 0.2), at.y = seq(0,0.4,by = 0.05),
xlab = 'Confidence score on ovu. estimates',ylab = '% of cycles' )
dev.off()## quartz_off_screen
## 27 HMM-state probabilities by cycle length
source("Scripts/_setup_Sympto.R")
load(file = paste0(IO$output_Rdata, 'states_probabilities_by_cycle_length_sympto.Rdata'), verbose = TRUE)## Loading objects:
## state_probsstate_probs_sympto = state_probs
g = ggplot(state_probs, aes(x = cycleday_from_ovu, y = prob, fill = state ))
g = g + geom_area(position = 'identity', alpha = 0.7) +
facet_grid(cycle_length_bin ~ ., scales = "free_y") +
scale_fill_manual(values = hmm_par$cycle.states.colors) +
xlab("Cycleday from estimated ovulation")+
ylab("Probabilities * Nb of cycles")
g = g + theme_hc() +
theme(strip.text.y = element_text(angle = 0),
strip.background = element_blank())
g
ggsave(filename = paste0(IO$panels,"S4_states_by_cycle_length_sympto.pdf"), plot = g, height = 5)## Saving 7 x 5 in imagesource("Scripts/_setup_Kindara.R")
load(file = paste0(IO$output_Rdata, 'states_probabilities_by_cycle_length_sympto.Rdata'), verbose = TRUE)## Loading objects:
## state_probsstate_probs_sympto = state_probs
load(file = paste0(IO$output_Rdata, 'states_probabilities_by_cycle_length_kindara.Rdata'), verbose = TRUE)## Loading objects:
## state_probsstate_probs_kindara = state_probs
g = ggplot(state_probs, aes(x = cycleday_from_ovu, y = prob, fill = state ))
g = g + geom_area(position = 'identity', alpha = 0.7) +
facet_grid(cycle_length_bin ~ ., scales = "free_y") +
scale_fill_manual(values = hmm_par$cycle.states.colors) +
xlab("Cycleday from estimated ovulation")+
ylab("Probabilities * Nb of cycles")
g = g + theme_hc() +
theme(strip.text.y = element_text(angle = 0),
strip.background = element_blank())
g
ggsave(filename = paste0(IO$panels,"S4_states_by_cycle_length_kindara.pdf"), plot = g, height = 5)## Saving 7 x 5 in imagestate_probs = combine.sum.state.probs(state_probs_kindara, state_probs_sympto)
g = ggplot(state_probs, aes(x = cycleday_from_ovu, y = prob, fill = state ))
g = g + geom_area(position = 'identity', alpha = 0.7) +
facet_grid(cycle_length_bin ~ ., scales = "fixed") +
scale_fill_manual(values = hmm_par$cycle.states.colors) +
xlab("Cycleday from estimated ovulation")+
ylab("Probabilities * Nb of cycles")
g = g + theme_hc() +
theme(strip.text.y = element_text(angle = 0),
strip.background = element_blank())
g
ggsave(filename = paste0(IO$panels,"4_states_by_cycle_length_combined.pdf"), plot = g, width = 7, height = 5)8 Figure S5 : comparison with previous studies
load(file = paste0(IO$output_Rdata,"all_studies_for_phases_comparison.Rdata") , verbose = TRUE)## Loading objects:
## all_studiesall = all_studies
g = ggplot(all, aes(col = study,x = mean, y = study)) + facet_grid( . ~ phase,scale = "free_x", space = "free_x")+
geom_point(shape = "|", size = 3, stroke = 6)+
guides(col = FALSE)+
scale_x_continuous(breaks = seq(8,46,by = 4))+
# geom_segment(aes(x = min, xend = max, yend = study), size = 0.5)+
geom_segment(aes(x = perc.2.5, xend = perc.97.5, yend = study), size = 2, alpha = 0.4)+
geom_segment(aes(x = perc.5, xend = perc.95, yend = study), size = 4, alpha = 0.4)+
geom_segment(aes(x = perc.10, xend = perc.90, yend = study), size = 5, alpha = 0.3)+
# geom_segment(aes(x = mean - SD, xend = mean + SD, yend = study), size = 0.5, alpha = 1)+
# geom_point(aes(x = mean - SD), shape = "|", size = 2.5)+ geom_point(aes(x = mean + SD), shape = "|", size = 2.5)+
xlab("duration")+ ylab("")+
theme_light()
g## Warning: Removed 6 rows containing missing values (geom_segment).## Warning: Removed 17 rows containing missing values (geom_segment).
## Warning: Removed 17 rows containing missing values (geom_segment).
ggsave(filename = paste0(IO$panels,"S5_cycle_phases_comparison_with_prev_studies.pdf"),plot = g, width = 6, height = 4)## Warning: Removed 6 rows containing missing values (geom_segment).
## Warning: Removed 17 rows containing missing values (geom_segment).
## Warning: Removed 17 rows containing missing values (geom_segment).g = ggplot(all, aes(col = study,x = range_95, y = study)) + facet_grid( . ~ phase,scale = "free_x", space = "free_x")+
geom_point(shape = "|", size = 3, stroke = 4)+
geom_point(aes(x = range_90), shape = "*", size = 3, stroke = 4)+
geom_point(aes(x = range_80), shape = 16, size = 1)+
geom_point(aes(x = 0), shape = 16, size = 0,alpha = 0)+
scale_x_continuous(breaks = seq(0,46,by = 4), expand = c(0,1))+
guides(col = FALSE)+
xlab("range")+ ylab("")+
theme_light()
g## Warning: Removed 6 rows containing missing values (geom_point).## Warning: Removed 17 rows containing missing values (geom_point).
## Warning: Removed 17 rows containing missing values (geom_point).
ggsave(filename = paste0(IO$panels,"S5_cycle_phases_range_comparison_with_prev_studies.pdf"),plot = g, width = 6, height = 4)## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 17 rows containing missing values (geom_point).
## Warning: Removed 17 rows containing missing values (geom_point).g = ggplot(all, aes(fill = study,y = range_95, x = study))+ coord_flip() + facet_grid( . ~ phase,scale = "free_x", space = "free_x")+
geom_bar(stat = "identity", alpha = 0.3, width = 0.5)+
geom_bar(aes(y = range_90),stat = "identity", alpha = 0.3, width = 0.75)+
geom_bar(aes(y = range_80),stat = "identity", alpha = 0.3, width = 1)+
scale_y_continuous(breaks = seq(0,46,by = 4))+
guides(fill = FALSE)+
theme_light()
g## Warning: Removed 6 rows containing missing values (position_stack).## Warning: Removed 17 rows containing missing values (position_stack).
## Warning: Removed 17 rows containing missing values (position_stack).
n_cycles = unique(all[,c("study","n.cycles")])
g = ggplot(n_cycles, aes(fill = study,y = n.cycles/1000, x = study)) + coord_flip()+
facet_grid( . ~ "# of cycles",scale = "free_x", space = "free_x")+
geom_bar(stat = "identity")+
guides(fill = FALSE)+
scale_y_continuous(breaks = seq(0,round(max(n_cycles$n.cycles/1000)/10)*10,len = 3))+
xlab("")+ ylab("(Thousands)")+
theme_light()
g
ggsave(filename = paste0(IO$panels,"S5_cycle_phases_comparison_with_prev_studies_n_cycles.pdf"),plot = g, width = 3, height = 4)9 Figure S6 : impact of noise and missing data on the estimations
9.1 Noise
load( file = paste0(IO$output_Rdata, "RES_1.Rdata"), verbose = TRUE)## Loading objects:
## RESxlim.max = 20
g.noise = ggplot(RES, aes(x = ovu.diff, fill = noise))
g.noise = g.noise + geom_histogram(bins = 100) +
facet_grid(noise.r ~.) +
scale_fill_gradient(low = "gray30", high = "tomato")+
guides(fill=FALSE)+
scale_y_continuous(breaks = seq(0,100,by =30))+
theme(strip.text.y = element_text(angle = 0))+
#ggtitle("difference in ovulation estimation") +
xlab("Diff. in ovu. est. (no noise - noise) [days]") + xlim(-xlim.max,xlim.max)
g.noise## Warning: Removed 3 rows containing non-finite values (stat_bin).
w = 15
ggsave(filename = paste0(IO$panels,"S6_noise_hist.pdf"), g.noise, width = w, height = w/2, units = "cm")## Warning: Removed 3 rows containing non-finite values (stat_bin).g = ggplot(RES, aes(x = dev.sym, fill = factor(noise)))
g + geom_histogram(bins = 50) + facet_grid(noise~.)
g.noise.scatter = ggplot(RES, aes(x = abs(ovu.diff), y = ovu.sd.no.noise, col = noise))
g.noise.scatter = g.noise.scatter +
geom_abline(slope = 1, intercept = 0, col = "gray")+
geom_point(size = 0.5) +
scale_color_gradient(low = "gray30", high = "tomato")+
facet_wrap(noise.r ~., ncol = 5)+
xlab("Abs. diff. in ovu. est. (no noise - noise) [days]")+
ylab("Uncertainty on ovulation est. (no noise)")+
guides(col=FALSE) + xlim(0,xlim.max)
g.noise.scatter## Warning: Removed 3 rows containing missing values (geom_point).
ggsave(paste0(IO$panels,"S6_noise_scatter.pdf"), g.noise.scatter, width = w, height = w/2, units = "cm")## Warning: Removed 3 rows containing missing values (geom_point).load( file = paste0(IO$output_Rdata, "RES_2.Rdata"), verbose = TRUE)## Loading objects:
## RES.cg.cons = ggplot(RES.c, aes(x = ovu.diff, fill = C))
g.cons = g.cons + geom_histogram(bins = 100) +
facet_grid(C ~.) +
scale_fill_gradient(low = "gray30", high = "slateblue1")+
guides(fill=FALSE)+
scale_y_continuous(breaks = seq(0,100,by =50))+
theme(strip.text.y = element_text(angle = 0))+
#ggtitle("difference in ovulation estimation") +
xlab("Diff. in ovu. est. (uniform - non-uniform) [days]") + xlim(-xlim.max, xlim.max)
g.cons
w = 15
ggsave(paste0(IO$panels,"S6_cons_hist.pdf"), g.cons, width = w, height = w/2, units = "cm")
g.cons.scatter = ggplot(RES.c, aes(x = abs(ovu.diff), y = ovu.sd.no.noise, col = C))
g.cons.scatter = g.cons.scatter +
geom_abline(slope = 1, intercept = 0, col = "gray")+
geom_point(size = 0.5) +
scale_color_gradient(low = "gray30", high = "slateblue1")+
facet_wrap(C ~., ncol = 5)+
xlab("Abs. diff. in ovu. est. (uniform - non-uniform) [days]")+
ylab("Uncertainty on ovulation est. (uniform)")+
guides(col=FALSE)+ xlim(0,xlim.max)
g.cons.scatter
ggsave(paste0(IO$panels,"S6_cons_scatter.pdf"), g.cons.scatter, width = w, height = w/2, units = "cm")load(file = paste0(IO$output_Rdata, "RES_3.Rdata"), verbose = TRUE)## Loading objects:
## RES.eg.em = ggplot(RES.e, aes(x = ovu.diff, fill = noise))
g.em = g.em + geom_histogram(bins = 100) +
facet_grid(noise.r ~.) +
scale_fill_gradient(low = "gray30", high = "green3")+
guides(fill=FALSE)+
scale_y_continuous(breaks = seq(0,100,by =30))+
theme(strip.text.y = element_text(angle = 0))+
#ggtitle("difference in ovulation estimation") +
xlab("Diff. in ovu. est. (no noise - noise) [days]") + xlim(-xlim.max, xlim.max)
g.em## Warning: Removed 3 rows containing non-finite values (stat_bin).
w = 15
ggsave(paste0(IO$panels,"S6_emission_hist.pdf"), g.em, width = w, height = w/2, units = "cm")## Warning: Removed 3 rows containing non-finite values (stat_bin).g.em.scatter = ggplot(RES.e, aes(x = abs(ovu.diff), y = ovu.sd.no.noise, col = noise.r))
g.em.scatter = g.em.scatter +
geom_abline(slope = 1, intercept = 0, col = "gray")+
geom_point(size = 0.5) +
scale_color_gradient(low = "gray30", high = "green3")+
facet_wrap(noise.r ~., ncol = 5)+
xlab("Abs. diff. in ovu. est. (no noise - noise) [days]")+
ylab("Uncertainty on ovulation est. (no noise)")+
guides(col=FALSE) + xlim(0,xlim.max)
g.em.scatter## Warning: Removed 3 rows containing missing values (geom_point).
w = 15
ggsave(paste0(IO$panels,"S6_emission_scatter.pdf"), g.em.scatter, width = w, height = w/2, units = "cm")## Warning: Removed 3 rows containing missing values (geom_point).g.em.sd = ggplot(RES.e, aes(x = ovu.sd, fill = noise))
g.em.sd = g.em.sd + geom_histogram(bins = 100) +
facet_grid(noise.r ~.) +
scale_fill_gradient(low = "gray30", high = "green3")+
guides(fill=FALSE)+
scale_y_continuous(breaks = seq(0,100,by =30))+
theme(strip.text.y = element_text(angle = 0))+
#ggtitle("difference in ovulation estimation") +
xlab("Ovu. est. uncertainty [days]")
g.em.sd
g.em.diff = ggplot(RES.e, aes(x = diff.em, fill = noise))
g.em.diff = g.em.diff + geom_histogram(bins = 100) +
facet_grid(noise.r ~.) +
scale_fill_gradient(low = "gray30", high = "green3")+
guides(fill=FALSE)+
scale_y_continuous(breaks = seq(0,100,by =30))+
theme(strip.text.y = element_text(angle = 0))+
#ggtitle("difference in ovulation estimation") +
xlab("Ovu. est. uncertainty [days]")
g.em.diff
9.2 missing data
load(file = paste0(IO$output_Rdata,"most_useful_sign_RES_l.Rdata"), verbose = TRUE)## Loading objects:
## RES_lg = ggplot(RES_l, aes(fill = percent_missing, col = percent_missing, y = value, x = percent_missing_f)) +
geom_hline(yintercept = 0, col = "gray70", size = 1)+
coord_flip() + geom_boxplot(outlier.size = 0.5, alpha = 0.5) +
facet_grid(missing_feature ~ metric, scale = "free")+
scale_fill_gradient(low = "gray20",high = "tomato") +
scale_color_gradient(low = "gray20",high = "tomato") +
guides(fill = FALSE, col = FALSE)+
ylab("comparison with complete profiles: difference [ovulation time] or log2(ratio) [sd and confidence]") + xlab("% of missing data")+
theme_light()
g
ggsave(filename = paste0(IO$panels,"S6_most_useful_sign_to_estimate_ovulation.pdf"), plot = g, width = 12, height = 3)