R Markdown
rm(list = ls())
library(lmerTest)
## 载入需要的程辑包：lme4
## 载入需要的程辑包：Matrix
##
## 载入程辑包：'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(splines)
library(lubridate)
library(sp)
library(gstat)
library(raster)
##
## 载入程辑包：'raster'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## The following object is masked from 'package:lme4':
##
## getData
library(dplyr)
setwd("D:/R/sperm quality and temperature")
data<-read.csv("D:/R/sperm quality and temperature/data_for_analysis.csv",fileEncoding = "GB18030")
data$daytemaver_cate_090<-cut(data$daytemaver_090, breaks = quantile(data$daytemaver_090,c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)),
labels = c("a","b","c","d","e","f","g","h","i","j"),right=FALSE,include.lowest=TRUE)
data$daytemaver_cate_09<-cut(data$daytemaver_09, breaks = quantile(data$daytemaver_09,c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)),
labels = c("a","b","c","d","e","f","g","h","i","j"),right=FALSE,include.lowest=TRUE)
data$daytemaver_cate_1014<-cut(data$daytemaver_1014, breaks = quantile(data$daytemaver_1014,c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)),
labels = c("a","b","c","d","e","f","g","h","i","j"),right=FALSE,include.lowest=TRUE)
data$daytemaver_cate_7090<-cut(data$daytemaver_7090, breaks = quantile(data$daytemaver_7090,c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)),
labels = c("a","b","c","d","e","f","g","h","i","j"),right=FALSE,include.lowest=TRUE)
##############################################################################
data$sumab<-data$sumab/sd(data$sumab)
data$sumabc<-data$sumabc/sd(data$sumabc)
data$density<-log(data$density)/sd(na.omit(log(data$density)))
data$semencount<-log(data$semencount)/sd(na.omit(log(data$semencount)))
data$totalmotilecount<-log(data$totalmotilecount)/sd(na.omit(log(data$totalmotilecount)))
data$progemotilecount<-log(data$progemotilecount)/sd(na.omit(log(data$progemotilecount)))
fit<-lmer(data=data,formula = sumab ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509+ gdp_aver + population_aver + urbanity+ ns(daytemaver_09,df=3))
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sumab ~ age.x + smokestatus + alcoholstatus + abstinenceday +
## bmi + occupation + education + season + green_nvdi + (1 |
## patientid) + hum_09 + pm2.509 + gdp_aver + population_aver +
## urbanity + ns(daytemaver_09, df = 3)
## Data: data
##
## REML criterion at convergence: 70652.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2172 -0.5880 -0.0195 0.5677 4.3002
##
## Random effects:
## Groups Name Variance Std.Dev.
## patientid (Intercept) 0.4408 0.6639
## Residual 0.5562 0.7458
## Number of obs: 26224, groups: patientid, 14063
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 2.537e+00 1.623e-01 2.172e+04 15.637
## age.x -1.160e-02 1.215e-03 2.288e+04 -9.552
## smokestatus 1.713e-01 1.618e-02 1.396e+04 10.592
## alcoholstatus -7.797e-02 9.132e-02 1.369e+04 -0.854
## abstinenceday -2.960e-02 5.327e-03 1.648e+04 -5.555
## bmi 6.213e-03 2.183e-03 1.387e+04 2.846
## occupationaaunemployees -1.098e-01 9.234e-02 1.430e+04 -1.189
## occupationagriculture -7.312e-02 1.037e-01 1.443e+04 -0.705
## occupationcivil_servant -1.285e-01 9.764e-02 1.406e+04 -1.316
## occupationmilitary_police -1.073e-01 1.084e-01 1.400e+04 -0.989
## occupationproduction_worker -9.919e-02 9.001e-02 1.426e+04 -1.102
## occupationsocial_production_service -6.898e-02 8.833e-02 1.427e+04 -0.781
## occupationtechnicist -1.389e-01 9.060e-02 1.418e+04 -1.533
## educationa初中以下 -2.399e-02 6.323e-02 1.386e+04 -0.379
## educationb高中专科 -1.786e-02 6.330e-02 1.382e+04 -0.282
## educationc本科以上 4.044e-02 6.466e-02 1.375e+04 0.626
## seasonbsummer -2.073e-02 2.472e-02 2.351e+04 -0.839
## seasoncautumn -4.032e-02 1.769e-02 2.524e+04 -2.279
## seasondwinter 4.925e-02 2.889e-02 2.419e+04 1.705
## green_nvdi -2.391e-02 5.460e-02 2.278e+04 -0.438
## hum_09 1.163e-04 7.613e-04 2.377e+04 0.153
## pm2.509 -1.109e-03 3.876e-04 2.473e+04 -2.862
## gdp_aver -2.411e-03 2.361e-03 1.371e+04 -1.021
## population_aver -8.749e-04 8.448e-04 1.344e+04 -1.036
## urbanity 5.637e-02 4.676e-02 1.376e+04 1.206
## ns(daytemaver_09, df = 3)1 9.865e-02 4.361e-02 2.358e+04 2.262
## ns(daytemaver_09, df = 3)2 -1.141e-01 1.285e-01 2.325e+04 -0.888
## ns(daytemaver_09, df = 3)3 -1.403e-01 5.204e-02 2.331e+04 -2.695
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## age.x < 2e-16 ***
## smokestatus < 2e-16 ***
## alcoholstatus 0.39327
## abstinenceday 2.81e-08 ***
## bmi 0.00443 **
## occupationaaunemployees 0.23453
## occupationagriculture 0.48056
## occupationcivil_servant 0.18820
## occupationmilitary_police 0.32252
## occupationproduction_worker 0.27048
## occupationsocial_production_service 0.43485
## occupationtechnicist 0.12523
## educationa初中以下 0.70439
## educationb高中专科 0.77789
## educationc本科以上 0.53163
## seasonbsummer 0.40163
## seasoncautumn 0.02268 *
## seasondwinter 0.08823 .
## green_nvdi 0.66142
## hum_09 0.87858
## pm2.509 0.00422 **
## gdp_aver 0.30717
## population_aver 0.30037
## urbanity 0.22802
## ns(daytemaver_09, df = 3)1 0.02370 *
## ns(daytemaver_09, df = 3)2 0.37458
## ns(daytemaver_09, df = 3)3 0.00704 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
fits<-lmer(data=data,formula = sumab ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi+ (1|patientid)
+ hum_09 + pm2.509+ gdp_aver + population_aver + urbanity + daytemaver_09)
comp1<-2*logLik(fit) - 2*logLik(fits)
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 1.548413e-05
pchisq_daytemaver_sumab_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_sumab_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_sumab_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.004213949 -0.1400247
## 2 27.15654 0.5303764 0.3144872 0.09248068 0.003454522 -0.1305504
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.002592973 -0.1380278
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.001421031 -0.1365823
## 5 26.67110 0.5512378 0.3142912 0.05921599 0.010200512 -0.1221453
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.006190162 -0.1424555
## yhci
## 1 0.1315968
## 2 0.1374594
## 3 0.1328418
## 4 0.1337402
## 5 0.1425463
## 6 0.1300752
daytemaver_sumab_09$miny<-daytemaver_sumab_09$y[daytemaver_sumab_09$y==min(daytemaver_sumab_09$y)]
daytemaver_sumab_09$minylci<-daytemaver_sumab_09$ylci[daytemaver_sumab_09$y==min(daytemaver_sumab_09$y)]
daytemaver_sumab_09$minyhci<-daytemaver_sumab_09$yhci[daytemaver_sumab_09$y==min(daytemaver_sumab_09$y)]
daytemaver_sumab_09$data.daytemaver_09[daytemaver_sumab_09$y==min(daytemaver_sumab_09$y)]
## [1] 34.86433
fit<-lmer(data=data,formula = sumabc ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509+ gdp_aver + population_aver + urbanity + ns(daytemaver_09,df=3))
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sumabc ~ age.x + smokestatus + alcoholstatus + abstinenceday +
## bmi + occupation + education + season + green_nvdi + (1 |
## patientid) + hum_09 + pm2.509 + gdp_aver + population_aver +
## urbanity + ns(daytemaver_09, df = 3)
## Data: data
##
## REML criterion at convergence: 70753.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8327 -0.5655 0.0561 0.6098 4.2807
##
## Random effects:
## Groups Name Variance Std.Dev.
## patientid (Intercept) 0.4031 0.6349
## Residual 0.5787 0.7607
## Number of obs: 26224, groups: patientid, 14063
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.256e+00 1.612e-01 2.190e+04 20.193
## age.x -1.040e-02 1.208e-03 2.263e+04 -8.607
## smokestatus 1.762e-01 1.593e-02 1.427e+04 11.056
## alcoholstatus -1.027e-01 8.991e-02 1.396e+04 -1.142
## abstinenceday -1.992e-02 5.263e-03 1.661e+04 -3.785
## bmi 6.402e-03 2.150e-03 1.418e+04 2.978
## occupationaaunemployees -1.665e-01 9.099e-02 1.467e+04 -1.830
## occupationagriculture -1.429e-01 1.022e-01 1.480e+04 -1.399
## occupationcivil_servant -1.583e-01 9.618e-02 1.441e+04 -1.646
## occupationmilitary_police -1.247e-01 1.068e-01 1.434e+04 -1.167
## occupationproduction_worker -1.255e-01 8.869e-02 1.463e+04 -1.415
## occupationsocial_production_service -9.406e-02 8.703e-02 1.464e+04 -1.081
## occupationtechnicist -1.738e-01 8.926e-02 1.454e+04 -1.948
## educationa初中以下 -2.503e-02 6.227e-02 1.413e+04 -0.402
## educationb高中专科 -8.854e-03 6.233e-02 1.408e+04 -0.142
## educationc本科以上 6.119e-02 6.366e-02 1.401e+04 0.961
## seasonbsummer -1.126e-02 2.491e-02 2.401e+04 -0.452
## seasoncautumn -4.982e-02 1.779e-02 2.552e+04 -2.801
## seasondwinter 4.554e-02 2.909e-02 2.462e+04 1.566
## green_nvdi -1.008e-02 5.430e-02 2.252e+04 -0.186
## hum_09 1.644e-04 7.670e-04 2.426e+04 0.214
## pm2.509 -2.641e-03 3.899e-04 2.512e+04 -6.772
## gdp_aver -2.155e-03 2.324e-03 1.399e+04 -0.927
## population_aver -4.665e-04 8.314e-04 1.370e+04 -0.561
## urbanity 3.230e-02 4.604e-02 1.404e+04 0.701
## ns(daytemaver_09, df = 3)1 1.691e-02 4.395e-02 2.408e+04 0.385
## ns(daytemaver_09, df = 3)2 -2.045e-01 1.296e-01 2.379e+04 -1.578
## ns(daytemaver_09, df = 3)3 -1.670e-01 5.246e-02 2.384e+04 -3.183
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## age.x < 2e-16 ***
## smokestatus < 2e-16 ***
## alcoholstatus 0.253295
## abstinenceday 0.000155 ***
## bmi 0.002910 **
## occupationaaunemployees 0.067250 .
## occupationagriculture 0.161811
## occupationcivil_servant 0.099805 .
## occupationmilitary_police 0.243153
## occupationproduction_worker 0.157043
## occupationsocial_production_service 0.279863
## occupationtechnicist 0.051473 .
## educationa初中以下 0.687734
## educationb高中专科 0.887051
## educationc本科以上 0.336434
## seasonbsummer 0.651391
## seasoncautumn 0.005096 **
## seasondwinter 0.117421
## green_nvdi 0.852762
## hum_09 0.830282
## pm2.509 1.29e-11 ***
## gdp_aver 0.353733
## population_aver 0.574742
## urbanity 0.483015
## ns(daytemaver_09, df = 3)1 0.700475
## ns(daytemaver_09, df = 3)2 0.114599
## ns(daytemaver_09, df = 3)3 0.001458 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
fits<-lmer(data=data,formula = sumabc ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi+ (1|patientid)
+ hum_09+ pm2.509 + gdp_aver + population_aver + urbanity+ daytemaver_09)
comp1<-2*logLik(fit) - 2*logLik(fits)
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 0.003792172
pchisq_daytemaver_sumabc_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
summary(pchisq_daytemaver_sumabc_09)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003792 0.003792 0.003792 0.003792 0.003792 0.003792
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_sumabc_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_sumabc_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.07745335 -0.2143580
## 2 27.15654 0.5303764 0.3144872 0.09248068 -0.07079080 -0.2058742
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.07605216 -0.2125777
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.07503690 -0.2112865
## 5 26.67110 0.5512378 0.3142912 0.05921599 -0.06484274 -0.1982530
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.07915707 -0.2165202
## yhci
## 1 0.05945133
## 2 0.06429264
## 3 0.06047338
## 4 0.06121272
## 5 0.06856751
## 6 0.05820612
daytemaver_sumabc_09$miny<-daytemaver_sumabc_09$y[daytemaver_sumabc_09$y==min(daytemaver_sumabc_09$y)]
daytemaver_sumabc_09$minylci<-daytemaver_sumabc_09$ylci[daytemaver_sumabc_09$y==min(daytemaver_sumabc_09$y)]
daytemaver_sumabc_09$minyhci<-daytemaver_sumabc_09$yhci[daytemaver_sumabc_09$y==min(daytemaver_sumabc_09$y)]
daytemaver_sumabc_09$data.daytemaver_09[daytemaver_sumabc_09$y==min(daytemaver_sumabc_09$y)]
## [1] 34.86433
fit<-lmer(data=data,formula = density ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity+ ns(daytemaver_09,df=3))
fits<-lmer(data=data,formula = density ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509+ gdp_aver + population_aver + urbanity + daytemaver_09)
comp1<-2*logLik(fit) - 2*logLik(fits)
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 7.105223e-06
pchisq_daytemaver_density_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_density_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_density_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 0.0027412112 -0.1268545
## 2 27.15654 0.5303764 0.3144872 0.09248068 0.0113415110 -0.1165332
## 3 27.55632 0.5103817 0.3155267 0.12068835 0.0045391977 -0.1246983
## 4 27.48101 0.5143356 0.3152724 0.11532073 0.0058453232 -0.1231314
## 5 26.67110 0.5512378 0.3142912 0.05921599 0.0191498430 -0.1071435
## 6 27.78178 0.4980376 0.3164460 0.13690500 0.0005618509 -0.1294671
## yhci
## 1 0.1323369
## 2 0.1392162
## 3 0.1337767
## 4 0.1348221
## 5 0.1454431
## 6 0.1305908
daytemaver_density_09$miny<-daytemaver_density_09$y[daytemaver_density_09$y==min(daytemaver_density_09$y)]
daytemaver_density_09$minylci<-daytemaver_density_09$ylci[daytemaver_density_09$y==min(daytemaver_density_09$y)]
daytemaver_density_09$minyhci<-daytemaver_density_09$yhci[daytemaver_density_09$y==min(daytemaver_density_09$y)]
daytemaver_density_09$data.daytemaver_09[daytemaver_density_09$y==min(daytemaver_density_09$y)]
## [1] 34.86433
fit<-lmer(data=data,formula = semencount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + ns(daytemaver_09,df=3))
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: semencount ~ age.x + smokestatus + alcoholstatus + abstinenceday +
## bmi + occupation + education + season + green_nvdi + (1 |
## patientid) + hum_09 + pm2.509 + gdp_aver + population_aver +
## urbanity + ns(daytemaver_09, df = 3)
## Data: data
##
## REML criterion at convergence: 69402.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2080 -0.4018 0.0897 0.5175 4.0019
##
## Random effects:
## Groups Name Variance Std.Dev.
## patientid (Intercept) 0.4828 0.6948
## Residual 0.5037 0.7097
## Number of obs: 26170, groups: patientid, 14054
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.365e+00 1.610e-01 2.196e+04 20.906
## age.x 8.829e-03 1.202e-03 2.367e+04 7.348
## smokestatus 5.061e-02 1.627e-02 1.446e+04 3.111
## alcoholstatus -8.738e-02 9.185e-02 1.424e+04 -0.951
## abstinenceday 5.841e-02 5.332e-03 1.719e+04 10.956
## bmi -4.316e-04 2.196e-03 1.438e+04 -0.197
## occupationaaunemployees -1.888e-01 9.277e-02 1.470e+04 -2.036
## occupationagriculture -7.962e-02 1.041e-01 1.483e+04 -0.765
## occupationcivil_servant -7.512e-02 9.815e-02 1.451e+04 -0.765
## occupationmilitary_police -2.152e-02 1.090e-01 1.447e+04 -0.197
## occupationproduction_worker -9.375e-02 9.044e-02 1.466e+04 -1.037
## occupationsocial_production_service -7.608e-02 8.875e-02 1.467e+04 -0.857
## occupationtechnicist -1.480e-01 9.104e-02 1.459e+04 -1.626
## educationa初中以下 -3.972e-02 6.364e-02 1.446e+04 -0.624
## educationb高中专科 4.357e-03 6.372e-02 1.441e+04 0.068
## educationc本科以上 9.723e-02 6.508e-02 1.435e+04 1.494
## seasonbsummer 6.475e-02 2.398e-02 2.299e+04 2.701
## seasoncautumn -2.141e-02 1.723e-02 2.483e+04 -1.243
## seasondwinter 6.879e-02 2.805e-02 2.371e+04 2.453
## green_nvdi 3.257e-03 5.404e-02 2.358e+04 0.060
## hum_09 -5.829e-04 7.385e-04 2.324e+04 -0.789
## pm2.509 -1.331e-03 3.768e-04 2.422e+04 -3.532
## gdp_aver 2.291e-03 2.375e-03 1.425e+04 0.965
## population_aver 1.508e-03 8.503e-04 1.402e+04 1.774
## urbanity -3.924e-02 4.704e-02 1.429e+04 -0.834
## ns(daytemaver_09, df = 3)1 6.140e-02 4.231e-02 2.307e+04 1.451
## ns(daytemaver_09, df = 3)2 -4.223e-02 1.246e-01 2.272e+04 -0.339
## ns(daytemaver_09, df = 3)3 -2.323e-01 5.056e-02 2.279e+04 -4.595
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## age.x 2.08e-13 ***
## smokestatus 0.001868 **
## alcoholstatus 0.341448
## abstinenceday < 2e-16 ***
## bmi 0.844162
## occupationaaunemployees 0.041817 *
## occupationagriculture 0.444470
## occupationcivil_servant 0.444105
## occupationmilitary_police 0.843558
## occupationproduction_worker 0.299955
## occupationsocial_production_service 0.391318
## occupationtechnicist 0.104072
## educationa初中以下 0.532582
## educationb高中专科 0.945480
## educationc本科以上 0.135235
## seasonbsummer 0.006928 **
## seasoncautumn 0.213888
## seasondwinter 0.014186 *
## green_nvdi 0.951939
## hum_09 0.429937
## pm2.509 0.000414 ***
## gdp_aver 0.334702
## population_aver 0.076095 .
## urbanity 0.404133
## ns(daytemaver_09, df = 3)1 0.146735
## ns(daytemaver_09, df = 3)2 0.734646
## ns(daytemaver_09, df = 3)3 4.36e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
fits<-lmer(data=data,formula = semencount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + daytemaver_09)
summary(fits)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: semencount ~ age.x + smokestatus + alcoholstatus + abstinenceday +
## bmi + occupation + education + season + green_nvdi + (1 |
## patientid) + hum_09 + pm2.509 + gdp_aver + population_aver +
## urbanity + daytemaver_09
## Data: data
##
## REML criterion at convergence: 69431
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2793 -0.4033 0.0907 0.5175 3.9794
##
## Random effects:
## Groups Name Variance Std.Dev.
## patientid (Intercept) 0.4825 0.6946
## Residual 0.5046 0.7104
## Number of obs: 26170, groups: patientid, 14054
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.471e+00 1.494e-01 1.988e+04 23.231
## age.x 8.913e-03 1.202e-03 2.366e+04 7.415
## smokestatus 5.011e-02 1.627e-02 1.446e+04 3.079
## alcoholstatus -9.106e-02 9.186e-02 1.424e+04 -0.991
## abstinenceday 5.829e-02 5.333e-03 1.719e+04 10.931
## bmi -3.304e-04 2.196e-03 1.438e+04 -0.150
## occupationaaunemployees -1.893e-01 9.279e-02 1.470e+04 -2.040
## occupationagriculture -7.957e-02 1.041e-01 1.483e+04 -0.764
## occupationcivil_servant -7.360e-02 9.817e-02 1.451e+04 -0.750
## occupationmilitary_police -2.015e-02 1.091e-01 1.447e+04 -0.185
## occupationproduction_worker -9.238e-02 9.046e-02 1.466e+04 -1.021
## occupationsocial_production_service -7.378e-02 8.876e-02 1.467e+04 -0.831
## occupationtechnicist -1.468e-01 9.106e-02 1.459e+04 -1.612
## educationa初中以下 -4.402e-02 6.365e-02 1.445e+04 -0.692
## educationb高中专科 5.642e-04 6.372e-02 1.441e+04 0.009
## educationc本科以上 9.320e-02 6.509e-02 1.435e+04 1.432
## seasonbsummer 7.916e-03 2.155e-02 2.304e+04 0.367
## seasoncautumn -2.317e-02 1.719e-02 2.480e+04 -1.348
## seasondwinter 1.098e-02 2.234e-02 2.375e+04 0.492
## green_nvdi 1.361e-02 5.368e-02 2.374e+04 0.254
## hum_09 -4.814e-04 7.287e-04 2.328e+04 -0.661
## pm2.509 -1.384e-03 3.722e-04 2.419e+04 -3.718
## gdp_aver 2.292e-03 2.375e-03 1.424e+04 0.965
## population_aver 1.496e-03 8.500e-04 1.401e+04 1.761
## urbanity -3.593e-02 4.704e-02 1.429e+04 -0.764
## daytemaver_09 -3.442e-03 1.438e-03 2.331e+04 -2.394
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## age.x 1.25e-13 ***
## smokestatus 0.002078 **
## alcoholstatus 0.321535
## abstinenceday < 2e-16 ***
## bmi 0.880399
## occupationaaunemployees 0.041391 *
## occupationagriculture 0.444832
## occupationcivil_servant 0.453448
## occupationmilitary_police 0.853424
## occupationproduction_worker 0.307160
## occupationsocial_production_service 0.405871
## occupationtechnicist 0.106899
## educationa初中以下 0.489168
## educationb高中专科 0.992936
## educationc本科以上 0.152222
## seasonbsummer 0.713340
## seasoncautumn 0.177784
## seasondwinter 0.622979
## green_nvdi 0.799843
## hum_09 0.508879
## pm2.509 0.000201 ***
## gdp_aver 0.334580
## population_aver 0.078332 .
## urbanity 0.444948
## daytemaver_09 0.016671 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 26 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
comp1<-2*logLik(fit) - 2*logLik(fits)
summary(comp1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.63 28.63 28.63 28.63 28.63 28.63
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 6.072744e-07
pchisq_daytemaver_semencount_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_semencount_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_semencount_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.012092721 -0.1437975
## 2 27.15654 0.5303764 0.3144872 0.09248068 -0.002203807 -0.1321514
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.010028534 -0.1413675
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.008528024 -0.1396008
## 5 26.67110 0.5512378 0.3142912 0.05921599 0.006813349 -0.1215195
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.014592696 -0.1467398
## yhci
## 1 0.1196120
## 2 0.1277438
## 3 0.1213105
## 4 0.1225448
## 5 0.1351462
## 6 0.1175544
daytemaver_semencount_09$miny<-daytemaver_semencount_09$y[daytemaver_semencount_09$y==min(daytemaver_semencount_09$y)]
daytemaver_semencount_09$minylci<-daytemaver_semencount_09$ylci[daytemaver_semencount_09$y==min(daytemaver_semencount_09$y)]
daytemaver_semencount_09$minyhci<-daytemaver_semencount_09$yhci[daytemaver_semencount_09$y==min(daytemaver_semencount_09$y)]
daytemaver_semencount_09$data.daytemaver_09[daytemaver_semencount_09$y==min(daytemaver_semencount_09$y)]
## [1] 34.86433
fit<-lmer(data=data,formula = totalmotilecount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + ns(daytemaver_09,df=3))
fits<-lmer(data=data,formula = totalmotilecount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + daytemaver_09)
comp1<-2*logLik(fit) - 2*logLik(fits)
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 2.471763e-07
pchisq_daytemaver_totalmotilecount_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_totalmotilecount_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_totalmotilecount_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.03009515 -0.1605124
## 2 27.15654 0.5303764 0.3144872 0.09248068 -0.01981176 -0.1484894
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.02794553 -0.1580007
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.02638390 -0.1561755
## 5 26.67110 0.5512378 0.3142912 0.05921599 -0.01047239 -0.1375514
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.03270056 -0.1635557
## yhci
## 1 0.10032211
## 2 0.10886592
## 3 0.10210963
## 4 0.10340773
## 5 0.11660664
## 6 0.09815458
daytemaver_totalmotilecount_09$miny<-daytemaver_totalmotilecount_09$y[daytemaver_totalmotilecount_09$y==min(daytemaver_totalmotilecount_09$y)]
daytemaver_totalmotilecount_09$minylci<-daytemaver_totalmotilecount_09$ylci[daytemaver_totalmotilecount_09$y==min(daytemaver_totalmotilecount_09$y)]
daytemaver_totalmotilecount_09$minyhci<-daytemaver_totalmotilecount_09$yhci[daytemaver_totalmotilecount_09$y==min(daytemaver_totalmotilecount_09$y)]
daytemaver_totalmotilecount_09$data.daytemaver_09[daytemaver_totalmotilecount_09$y==min(daytemaver_totalmotilecount_09$y)]
## [1] 34.86433
fit<-lmer(data=data,formula = progemotilecount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + ns(daytemaver_09,df=3))
fits<-lmer(data=data,formula = progemotilecount ~ age.x + smokestatus + alcoholstatus + abstinenceday + bmi + occupation + education + season + green_nvdi + (1|patientid)
+ hum_09 + pm2.509 + gdp_aver + population_aver + urbanity + daytemaver_09)
comp1<-2*logLik(fit) - 2*logLik(fits)
pchisq(as.numeric(comp1), df=2, lower.tail=F)
## [1] 9.440661e-09
pchisq_daytemaver_progemotilecount_09<-pchisq(as.numeric(comp1), df=2, lower.tail=F)
a<-data.frame(ns(data$daytemaver_09,df=3))
y<-(summary(fit)$coefficients[26,1])*(a$X1)+ (summary(fit)$coefficients[27,1])*(a$X2)+ (summary(fit)$coefficients[28,1])*(a$X3)
ylci<-(summary(fit)$coefficients[26,1]-1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]-1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]-1.96*summary(fit)$coefficients[28,2])*(a$X3)
yhci<-(summary(fit)$coefficients[26,1]+1.96*summary(fit)$coefficients[26,2])*(a$X1)+ (summary(fit)$coefficients[27,1]+1.96*summary(fit)$coefficients[27,2])*(a$X2)+ (summary(fit)$coefficients[28,1]+1.96*summary(fit)$coefficients[28,2])*(a$X3)
daytemaver_progemotilecount_09<-data.frame(data$daytemaver_09,a,y,ylci,yhci)
head(daytemaver_progemotilecount_09)
## data.daytemaver_09 X1 X2 X3 y ylci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.004889477 -0.1347204
## 2 27.15654 0.5303764 0.3144872 0.09248068 0.005784644 -0.1223147
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.002653018 -0.1321235
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.001029915 -0.1302381
## 5 26.67110 0.5512378 0.3142912 0.05921599 0.015416052 -0.1110920
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.007603469 -0.1378703
## yhci
## 1 0.1249415
## 2 0.1338840
## 3 0.1268175
## 4 0.1281783
## 5 0.1419241
## 6 0.1226634
daytemaver_progemotilecount_09$miny<-daytemaver_progemotilecount_09$y[daytemaver_progemotilecount_09$y==min(daytemaver_progemotilecount_09$y)]
daytemaver_progemotilecount_09$minylci<-daytemaver_progemotilecount_09$ylci[daytemaver_progemotilecount_09$y==min(daytemaver_progemotilecount_09$y)]
daytemaver_progemotilecount_09$minyhci<-daytemaver_progemotilecount_09$yhci[daytemaver_progemotilecount_09$y==min(daytemaver_progemotilecount_09$y)]
daytemaver_progemotilecount_09$data.daytemaver_09[daytemaver_progemotilecount_09$y==min(daytemaver_progemotilecount_09$y)]
## [1] 34.86433
daytemaver_sumab_09$outcome<-"asumab"
daytemaver_sumabc_09$outcome<-"bsumabc"
daytemaver_density_09$outcome<-"cdensity"
daytemaver_semencount_09$outcome<-"dsemencount"
daytemaver_totalmotilecount_09$outcome<-"etotalmotilecount"
daytemaver_progemotilecount_09$outcome<-"fprogemotilecount"
##############################################################################
names(daytemaver_sumab_09)[1]<-"daytemaver"
names(daytemaver_sumabc_09)[1]<-"daytemaver"
names(daytemaver_density_09)[1]<-"daytemaver"
names(daytemaver_semencount_09)[1]<-"daytemaver"
names(daytemaver_totalmotilecount_09)[1]<-"daytemaver"
names(daytemaver_progemotilecount_09)[1]<-"daytemaver"
dataplot<-rbind(daytemaver_sumab_09,daytemaver_sumabc_09,daytemaver_density_09,daytemaver_semencount_09,daytemaver_totalmotilecount_09,daytemaver_progemotilecount_09)
head(dataplot,10)
## daytemaver X1 X2 X3 y ylci yhci
## 1 27.65894 0.5048565 0.3159160 0.12804251 -0.0042139490 -0.1400247 0.1315968
## 2 27.15654 0.5303764 0.3144872 0.09248068 0.0034545224 -0.1305504 0.1374594
## 3 27.55632 0.5103817 0.3155267 0.12068835 -0.0025929731 -0.1380278 0.1328418
## 4 27.48101 0.5143356 0.3152724 0.11532073 -0.0014210314 -0.1365823 0.1337402
## 5 26.67110 0.5512378 0.3142912 0.05921599 0.0102005122 -0.1221453 0.1425463
## 6 27.78178 0.4980376 0.3164460 0.13690500 -0.0061901624 -0.1424555 0.1300752
## 7 27.41477 0.5177426 0.3150709 0.11061963 -0.0004025654 -0.1353247 0.1345196
## 8 27.35818 0.5205997 0.3149153 0.10661924 0.0004581275 -0.1342611 0.1351773
## 9 27.07366 0.5342087 0.3143692 0.08672352 0.0046535324 -0.1290620 0.1383691
## 10 26.79657 0.5462169 0.3142261 0.06770651 0.0085217631 -0.1242445 0.1412880
## miny minylci minyhci outcome
## 1 -0.1637932 -0.3319094 0.004323059 asumab
## 2 -0.1637932 -0.3319094 0.004323059 asumab
## 3 -0.1637932 -0.3319094 0.004323059 asumab
## 4 -0.1637932 -0.3319094 0.004323059 asumab
## 5 -0.1637932 -0.3319094 0.004323059 asumab
## 6 -0.1637932 -0.3319094 0.004323059 asumab
## 7 -0.1637932 -0.3319094 0.004323059 asumab
## 8 -0.1637932 -0.3319094 0.004323059 asumab
## 9 -0.1637932 -0.3319094 0.004323059 asumab
## 10 -0.1637932 -0.3319094 0.004323059 asumab
pvalue<-c(pchisq_daytemaver_sumab_09,pchisq_daytemaver_sumabc_09,pchisq_daytemaver_density_09,pchisq_daytemaver_semencount_09,pchisq_daytemaver_totalmotilecount_09,pchisq_daytemaver_progemotilecount_09)
pdate<-c("a09","a09","a09","a09","a09","a09")
ptype<-c("asumab","bsumabc","cdensity","dsemencount","etotalmotilecount","fprogemotilecount")
ptext<-data.frame(pvalue=pvalue,outcome=ptype)
ptext$pvalue<-format(ptext$pvalue,scientific=T,digits=3)
ptext$label<-paste("P non-linear=",ptext$pvalue)
ptext
## pvalue outcome label
## 1 1.55e-05 asumab P non-linear= 1.55e-05
## 2 3.79e-03 bsumabc P non-linear= 3.79e-03
## 3 7.11e-06 cdensity P non-linear= 7.11e-06
## 4 6.07e-07 dsemencount P non-linear= 6.07e-07
## 5 2.47e-07 etotalmotilecount P non-linear= 2.47e-07
## 6 9.44e-09 fprogemotilecount P non-linear= 9.44e-09
outcome<-as_labeller(c("asumab"="Progressive motility","bsumabc"="Total motility","cdensity"="Sperm concentration",
"dsemencount"="Sperm count","etotalmotilecount"="Total motile sperm count","fprogemotilecount"="Progressive motile sperm count"))
plot_daytemaver_non_linear<- ggplot(dataplot, aes(daytemaver, y-miny)) +
geom_line(aes(daytemaver, y-miny)) +
geom_ribbon(aes(ymin = ylci-minylci, ymax = yhci-minyhci),fill='grey',alpha=0.5) +
geom_hline(aes(yintercept=0),color="#00BFC4",linetype="dashed") +
geom_line(aes(daytemaver,ylci-minylci),color="#CC6600") +
geom_line(aes(daytemaver,yhci-minyhci),color="blue") +
labs(x = "Daytemaver exposure",y = "Regression coefficient") +
facet_grid(cols =vars(outcome),scales = "free",labeller = labeller(outcome=outcome)) +
geom_text(data=ptext,mapping=aes(x=Inf, y=Inf,label=label),size=3, hjust="right",vjust="top",color="#9933FF")
tiff(filename="D:/R/sperm quality and temperature/daytemaver_non_linear/daytemaver_non_linear09.tiff",width=14.4,height=3.29,unit="in",res=900,compression="lzw")
plot_daytemaver_non_linear
dev.off()
## png
## 2