lm-antibodies
ztc
2022/1/25
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
data1 <- read.csv("D:/R/GWAS_ALL1-1_m.csv",head=T)
data1$male <- as.factor(as.character(data1$male))
data1$nation <- as.factor(as.character(data1$nation) )
data1$hypertension <- as.factor(as.character(data1$hypertension))
data1$diabetes <- as.factor(as.character(data1$diabetes))
data1$bands <- as.factor(as.character(data1$bands))
#data1$male <- as.factor(data1$male)
library(lme4)## 载入需要的程辑包：Matrix# compute p value by anova----
fit_NAb <- lmer(log(NAb) ~ nation + BMI + bands + dose_inter_days +
agec + inter_S_days + male +
diabetes + hypertension + NOD + (1 | num),
data = data1
)
fit_NAb_inter_S_days <- lmer(log(NAb) ~ nation + BMI + bands + dose_inter_days +
(agec + inter_S_days + male)^2 +
diabetes + hypertension + NOD + (1 | num),
data = data1
)
anova(fit_NAb, fit_NAb_inter_S_days)## refitting model(s) with ML (instead of REML)## Data: data1
## Models:
## fit_NAb: log(NAb) ~ nation + BMI + bands + dose_inter_days + agec + inter_S_days + male + diabetes + hypertension + NOD + (1 | num)
## fit_NAb_inter_S_days: log(NAb) ~ nation + BMI + bands + dose_inter_days + (agec + inter_S_days + male)^2 + diabetes + hypertension + NOD + (1 | num)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## fit_NAb 17 14525 14636 -7245.4 14491
## fit_NAb_inter_S_days 24 14515 14672 -7233.4 14467 24.093 7 0.001097 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# check interaction term----
fit_NAb <- lmerTest::lmer(log(NAb) ~nation + BMI + bands + dose_inter_days +
(agec + inter_S_days + male)^2 +
diabetes + hypertension + NOD + (1 | num),
data = data1
)
summary(fit_NAb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log(NAb) ~ nation + BMI + bands + dose_inter_days + (agec + inter_S_days +
## male)^2 + diabetes + hypertension + NOD + (1 | num)
## Data: data1
##
## REML criterion at convergence: 14615.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.5887 -0.4593 -0.0239 0.4570 5.1671
##
## Random effects:
## Groups Name Variance Std.Dev.
## num (Intercept) 0.2654 0.5151
## Residual 0.6668 0.8166
## Number of obs: 5264, groups: num, 4339
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.127e+00 1.268e-01 4.633e+03 16.772 < 2e-16 ***
## nation1 1.599e-02 3.603e-02 4.080e+03 0.444 0.65716
## BMI -1.942e-02 4.925e-03 4.245e+03 -3.943 8.18e-05 ***
## bands1 4.292e-01 3.557e-02 4.001e+03 12.066 < 2e-16 ***
## bands2 6.327e-01 3.452e-02 3.945e+03 18.331 < 2e-16 ***
## dose_inter_days 1.240e-02 1.400e-03 4.512e+03 8.861 < 2e-16 ***
## agec28-37 -6.234e-02 5.902e-02 5.019e+03 -1.056 0.29092
## agec38-47 -1.514e-01 6.882e-02 5.148e+03 -2.200 0.02784 *
## agec47+ -3.210e-01 7.333e-02 5.165e+03 -4.378 1.22e-05 ***
## inter_S_days -9.490e-03 4.979e-04 2.984e+03 -19.061 < 2e-16 ***
## male1 -1.044e-01 6.723e-02 5.242e+03 -1.553 0.12038
## diabetes1 -4.369e-01 1.622e-01 4.866e+03 -2.694 0.00709 **
## hypertension1 -9.241e-02 9.324e-02 4.670e+03 -0.991 0.32170
## NOD1 8.351e-02 9.286e-02 4.693e+03 0.899 0.36855
## NOD1+ -1.248e-01 1.846e-01 4.726e+03 -0.676 0.49895
## agec28-37:inter_S_days 1.429e-04 6.329e-04 2.853e+03 0.226 0.82133
## agec38-47:inter_S_days 4.126e-04 7.356e-04 2.939e+03 0.561 0.57494
## agec47+:inter_S_days 2.467e-03 8.151e-04 2.922e+03 3.027 0.00249 **
## agec28-37:male1 3.235e-02 8.155e-02 4.262e+03 0.397 0.69163
## agec38-47:male1 -1.410e-01 8.976e-02 4.337e+03 -1.571 0.11620
## agec47+:male1 -2.314e-01 9.719e-02 4.342e+03 -2.381 0.01732 *
## inter_S_days:male1 -1.462e-03 6.006e-04 3.308e+03 -2.434 0.01497 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itconfint(fit_NAb, oldNames=FALSE)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|num 0.4505825701 0.5675424479
## sigma 0.7806924322 0.8520583335
## (Intercept) 1.8786215288 2.3746678082
## nation1 -0.0544797602 0.0864701017
## BMI -0.0290514226 -0.0097857828
## bands1 0.3596199125 0.4987864354
## bands2 0.5651867598 0.7002431025
## dose_inter_days 0.0096666220 0.0151430768
## agec28-37 -0.1778967812 0.0531003505
## agec38-47 -0.2860339077 -0.0167817588
## agec47+ -0.4644991841 -0.1775826774
## inter_S_days -0.0104643969 -0.0085156186
## male1 -0.2359299450 0.0271558065
## diabetes1 -0.7541977263 -0.1196789476
## hypertension1 -0.2748402010 0.0899614471
## NOD1 -0.0981449083 0.2652759957
## NOD1+ -0.4859403294 0.2364556333
## agec28-37:inter_S_days -0.0010953700 0.0013865031
## agec38-47:inter_S_days -0.0010265142 0.0018530750
## agec47+:inter_S_days 0.0008727487 0.0040641902
## agec28-37:male1 -0.1271658406 0.1918758926
## agec38-47:male1 -0.3166481721 0.0345322393
## agec47+:male1 -0.4215184487 -0.0412976876
## inter_S_days:male1 -0.0026392937 -0.0002874007plot(fit_NAb)
which( residuals(fit_NAb)< -4 )## 1116 3472 4438
## 1114 3463 4424which( residuals(fit_NAb)> 4 )## 2112
## 2108hist(residuals(fit_NAb))
# qqnorm(y = residuals(fit_NAb))
# qqline(y = residuals(fit_NAb))
#
# fit_NAb_2 <- lmerTest::lmer(log(NAb) ~nation + BMI + bands + dose_inter_days +
# (agec + inter_S_days + male)^2 +
# diabetes + hypertension + NOD + (1 | num),
# data = data1[!(residuals(fit_NAb) < -4 | residuals(fit_NAb)> 4),],
#
# )
# plot(fit_NAb_2)
# which( residuals(fit_NAb_2)< -4 )
# which( residuals(fit_NAb_2)> 4 )
## calculate effects----
library(effects)## 载入需要的程辑包：carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.effect_fit_NAb <- allEffects(fit_NAb, confidence.level=.95)
effect_fit_NAb$`inter_S_days:male`##
## inter_S_days*male effect
## male
## inter_S_days 0 1
## 1 2.3673278 2.2122161
## 52 1.9074155 1.6777331
## 100 1.4745568 1.1746903
## 150 1.0236624 0.6506873
## 200 0.5727679 0.1266844effect_fit_NAb[["agec:male"]][["fit"]][5:8]/effect_fit_NAb[["agec:male"]][["fit"]][1:4]## [1] 0.8946851 0.9094099 0.8062325 0.7486758effect_fit_NAb[["agec:male"]][["upper"]][5:8]/effect_fit_NAb[["agec:male"]][["upper"]][1:4]## [1] 0.9179000 0.9386189 0.8379834 0.7854585effect_fit_NAb[["agec:male"]][["lower"]][5:8]/effect_fit_NAb[["agec:male"]][["lower"]][1:4]## [1] 0.8700750 0.8782736 0.7715118 0.7077326effect_fit_NAb$`agec:male`$confidence.level## [1] 0.95# marginal effect plot----
class(effect_fit_NAb)## [1] "efflist"plot(effect_fit_NAb)
plot(effect_fit_NAb,1)
plot(effect_fit_NAb,9)
plot(effect_fit_NAb,8)
plot(effect_fit_NAb,cex = 0.5)
# tran month to factor----
fit_NAb_mf <- lmerTest::lmer(log(NAb) ~nation + BMI + bands + dose_inter_days +
(agec + factor(inter_S_month) + male)^2 +
diabetes + hypertension + NOD + (1 | num),
data = data1
)
summary(fit_NAb_mf)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log(NAb) ~ nation + BMI + bands + dose_inter_days + (agec + factor(inter_S_month) +
## male)^2 + diabetes + hypertension + NOD + (1 | num)
## Data: data1
##
## REML criterion at convergence: 14478.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4719 -0.4588 -0.0152 0.4618 4.7897
##
## Random effects:
## Groups Name Variance Std.Dev.
## num (Intercept) 0.2818 0.5309
## Residual 0.6330 0.7956
## Number of obs: 5264, groups: num, 4339
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.875e+00 1.264e-01 4.509e+03 14.835
## nation1 1.548e-02 3.584e-02 4.037e+03 0.432
## BMI -1.947e-02 4.896e-03 4.199e+03 -3.978
## bands1 5.105e-01 3.681e-02 4.454e+03 13.870
## bands2 7.198e-01 3.606e-02 4.468e+03 19.962
## dose_inter_days 1.689e-02 1.542e-03 4.799e+03 10.952
## agec28-37 -7.588e-02 5.686e-02 5.208e+03 -1.335
## agec38-47 -2.029e-01 6.499e-02 5.164e+03 -3.122
## agec47+ -3.379e-01 6.904e-02 5.117e+03 -4.894
## factor(inter_S_month)2 -4.933e-01 9.005e-02 5.214e+03 -5.478
## factor(inter_S_month)3 -9.099e-01 6.763e-02 5.077e+03 -13.454
## factor(inter_S_month)4 -1.066e+00 5.926e-02 3.966e+03 -17.991
## factor(inter_S_month)5 -1.090e+00 9.067e-02 4.677e+03 -12.018
## factor(inter_S_month)6 -9.977e-01 1.360e-01 3.068e+03 -7.334
## factor(inter_S_month)7 -1.393e+00 4.607e-01 3.614e+03 -3.023
## male1 -1.796e-01 6.828e-02 5.119e+03 -2.631
## diabetes1 -3.926e-01 1.612e-01 4.814e+03 -2.435
## hypertension1 -3.649e-02 9.294e-02 4.622e+03 -0.393
## NOD1 7.738e-02 9.253e-02 4.643e+03 0.836
## NOD1+ -2.035e-01 1.835e-01 4.676e+03 -1.109
## agec28-37:factor(inter_S_month)2 -7.490e-02 1.348e-01 5.202e+03 -0.556
## agec38-47:factor(inter_S_month)2 1.826e-01 1.591e-01 5.213e+03 1.147
## agec47+:factor(inter_S_month)2 4.390e-02 1.544e-01 5.214e+03 0.284
## agec28-37:factor(inter_S_month)3 -4.392e-02 9.998e-02 4.943e+03 -0.439
## agec38-47:factor(inter_S_month)3 8.036e-02 1.344e-01 5.024e+03 0.598
## agec47+:factor(inter_S_month)3 -7.779e-02 1.771e-01 5.200e+03 -0.439
## agec28-37:factor(inter_S_month)4 -7.319e-02 8.380e-02 3.796e+03 -0.873
## agec38-47:factor(inter_S_month)4 -2.311e-02 9.580e-02 3.689e+03 -0.241
## agec47+:factor(inter_S_month)4 1.178e-01 1.118e-01 3.761e+03 1.053
## agec28-37:factor(inter_S_month)5 -6.004e-02 1.278e-01 4.431e+03 -0.470
## agec38-47:factor(inter_S_month)5 -4.060e-02 1.416e-01 4.058e+03 -0.287
## agec47+:factor(inter_S_month)5 5.680e-01 1.604e-01 3.863e+03 3.542
## agec28-37:factor(inter_S_month)6 -9.224e-02 1.615e-01 3.055e+03 -0.571
## agec38-47:factor(inter_S_month)6 9.848e-02 1.907e-01 3.331e+03 0.516
## agec47+:factor(inter_S_month)6 1.427e-01 2.007e-01 3.390e+03 0.711
## agec28-37:factor(inter_S_month)7 -1.463e-01 5.421e-01 4.083e+03 -0.270
## agec38-47:factor(inter_S_month)7 -6.084e-01 6.554e-01 4.738e+03 -0.928
## agec47+:factor(inter_S_month)7 -2.191e-02 1.122e+00 3.266e+03 -0.020
## agec28-37:male1 5.358e-02 8.339e-02 4.334e+03 0.643
## agec38-47:male1 -1.019e-01 9.110e-02 4.337e+03 -1.118
## agec47+:male1 -1.811e-01 9.896e-02 4.362e+03 -1.830
## factor(inter_S_month)2:male1 -4.242e-02 1.274e-01 5.205e+03 -0.333
## factor(inter_S_month)3:male1 -7.430e-03 9.641e-02 5.167e+03 -0.077
## factor(inter_S_month)4:male1 -1.953e-04 8.086e-02 4.270e+03 -0.002
## factor(inter_S_month)5:male1 -3.738e-01 1.228e-01 4.351e+03 -3.044
## factor(inter_S_month)6:male1 -1.926e-01 1.558e-01 3.571e+03 -1.237
## factor(inter_S_month)7:male1 -3.285e-02 4.657e-01 4.656e+03 -0.071
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## nation1 0.665778
## BMI 7.07e-05 ***
## bands1 < 2e-16 ***
## bands2 < 2e-16 ***
## dose_inter_days < 2e-16 ***
## agec28-37 0.182072
## agec38-47 0.001806 **
## agec47+ 1.02e-06 ***
## factor(inter_S_month)2 4.50e-08 ***
## factor(inter_S_month)3 < 2e-16 ***
## factor(inter_S_month)4 < 2e-16 ***
## factor(inter_S_month)5 < 2e-16 ***
## factor(inter_S_month)6 2.83e-13 ***
## factor(inter_S_month)7 0.002522 **
## male1 0.008544 **
## diabetes1 0.014909 *
## hypertension1 0.694634
## NOD1 0.403077
## NOD1+ 0.267353
## agec28-37:factor(inter_S_month)2 0.578415
## agec38-47:factor(inter_S_month)2 0.251343
## agec47+:factor(inter_S_month)2 0.776202
## agec28-37:factor(inter_S_month)3 0.660465
## agec38-47:factor(inter_S_month)3 0.549941
## agec47+:factor(inter_S_month)3 0.660499
## agec28-37:factor(inter_S_month)4 0.382513
## agec38-47:factor(inter_S_month)4 0.809387
## agec47+:factor(inter_S_month)4 0.292487
## agec28-37:factor(inter_S_month)5 0.638421
## agec38-47:factor(inter_S_month)5 0.774347
## agec47+:factor(inter_S_month)5 0.000402 ***
## agec28-37:factor(inter_S_month)6 0.567839
## agec38-47:factor(inter_S_month)6 0.605669
## agec47+:factor(inter_S_month)6 0.477281
## agec28-37:factor(inter_S_month)7 0.787269
## agec38-47:factor(inter_S_month)7 0.353306
## agec47+:factor(inter_S_month)7 0.984424
## agec28-37:male1 0.520539
## agec38-47:male1 0.263574
## agec47+:male1 0.067256 .
## factor(inter_S_month)2:male1 0.739137
## factor(inter_S_month)3:male1 0.938574
## factor(inter_S_month)4:male1 0.998073
## factor(inter_S_month)5:male1 0.002346 **
## factor(inter_S_month)6:male1 0.216323
## factor(inter_S_month)7:male1 0.943758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 47 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# use inter_S_month to fit lmer----
fit_NAb_m <- lmerTest::lmer(log(NAb) ~nation + BMI + bands + dose_inter_days +
(agec + inter_S_month + male)^2 +
diabetes + hypertension + NOD + (1 | num),
data = data1
)
summary(fit_NAb_m)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log(NAb) ~ nation + BMI + bands + dose_inter_days + (agec + inter_S_month +
## male)^2 + diabetes + hypertension + NOD + (1 | num)
## Data: data1
##
## REML criterion at convergence: 14595.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4947 -0.4628 -0.0249 0.4595 5.1901
##
## Random effects:
## Groups Name Variance Std.Dev.
## num (Intercept) 0.2806 0.5297
## Residual 0.6557 0.8097
## Number of obs: 5264, groups: num, 4339
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.279e+00 1.297e-01 4.749e+03 17.573 < 2e-16 ***
## nation1 1.728e-02 3.617e-02 4.073e+03 0.478 0.632915
## BMI -1.951e-02 4.943e-03 4.235e+03 -3.948 8.01e-05 ***
## bands1 4.194e-01 3.574e-02 4.001e+03 11.736 < 2e-16 ***
## bands2 6.151e-01 3.471e-02 3.954e+03 17.720 < 2e-16 ***
## dose_inter_days 1.223e-02 1.404e-03 4.497e+03 8.712 < 2e-16 ***
## agec28-37 -6.023e-02 6.620e-02 4.703e+03 -0.910 0.363032
## agec38-47 -1.521e-01 7.678e-02 4.864e+03 -1.980 0.047708 *
## agec47+ -3.721e-01 8.111e-02 4.908e+03 -4.588 4.58e-06 ***
## inter_S_month -2.866e-01 1.515e-02 3.078e+03 -18.926 < 2e-16 ***
## male1 -7.604e-02 7.268e-02 5.184e+03 -1.046 0.295530
## diabetes1 -4.235e-01 1.626e-01 4.845e+03 -2.604 0.009237 **
## hypertension1 -9.408e-02 9.352e-02 4.653e+03 -1.006 0.314502
## NOD1 7.619e-02 9.315e-02 4.676e+03 0.818 0.413443
## NOD1+ -1.341e-01 1.852e-01 4.708e+03 -0.724 0.469129
## agec28-37:inter_S_month 4.820e-03 1.920e-02 2.924e+03 0.251 0.801822
## agec38-47:inter_S_month 1.240e-02 2.227e-02 2.960e+03 0.557 0.577755
## agec47+:inter_S_month 8.089e-02 2.454e-02 2.976e+03 3.297 0.000988 ***
## agec28-37:male1 3.392e-02 8.181e-02 4.250e+03 0.415 0.678451
## agec38-47:male1 -1.402e-01 9.004e-02 4.324e+03 -1.557 0.119587
## agec47+:male1 -2.230e-01 9.752e-02 4.330e+03 -2.287 0.022270 *
## inter_S_month:male1 -4.712e-02 1.825e-02 3.316e+03 -2.583 0.009850 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## calculate effects----
effect_fit_NAb_m <- allEffects(fit_NAb_m, confidence.level=.95)
#age vs log NAb
# ggplot(data1,aes(x = inter_S_month, y = log(NAb),color = agec))+
# geom_point()+
# geom_jitter()
# calculate mean by month and agec----
library(dplyr)##
## 载入程辑包：'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionNAb_summ_by_month <- data1 %>%
group_by(inter_S_month,agec) %>%
summarise(mean_NAb = mean(log(NAb)),
se_NAb = sd(log(NAb))/sqrt(n()))# %>% ## `summarise()` has grouped output by 'inter_S_month'. You can override using the `.groups` argument.#mutate(inter_S_days = rep(15+(0:6)*30,each = 4))
NAb_summ_by_month$inter_S_days <- rep(15+(0:6)*30,each = 4)
# ggplot on mean by month and agec----
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v stringr 1.4.0
## v tidyr 1.1.4 v forcats 0.5.1
## v readr 2.1.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()ggplot(data1,aes(x = inter_S_days, y = log(NAb),color = agec))+
geom_point(alpha = 0.2)+
geom_errorbar(data = NAb_summ_by_month,
aes( y = mean_NAb,
ymin = mean_NAb - se_NAb,
ymax = mean_NAb + se_NAb,
),width=1,size = 1)+
geom_line(data = NAb_summ_by_month,
aes( y = mean_NAb),size = 1)+
theme_classic()