DNase_eda
Дмитрий Пасько
24 04 2020
Подключение библиотек + написание дополнительных функций:
library(tidyverse)## -- Attaching packages ---------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.1
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## Warning: package 'stringr' was built under R version 4.0.0## -- Conflicts ------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(magrittr)##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some# функция для оценки качества модели
mysummary = function(mdl) {
cat("-----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:\n")
cat("-----> Ошибка = " ,mdl$residuals^2 %>% mean(),'\n')
cat("\n")
gvlma::gvlma(mdl) %>% summary()
cat("\n")
cat("-----> БАЗОВЫЕ ГРАФИКИ:\n")
cat("\n")
par(mfrow = c(2, 2))
plot(mdl)
par(mfrow = c(1, 1))
cat("\n")
cat("-----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ\n")
cat("\n")
shapiro.test(mdl$residuals) %>% print()
cat("\n")
qqPlot(mdl, main = "Q-Q plot")
}Берётся такой набор данных:
df=DNaseОсновные статистики:
# основные статистики
psych::describe(df) ## vars n mean sd median trimmed mad min max range skew kurtosis
## Run* 1 176 6.00 3.17 6.00 6.00 4.45 1.00 11.0 10.00 0.00 -1.24
## conc 2 176 3.11 4.06 1.17 2.35 1.56 0.05 12.5 12.45 1.44 0.74
## density 3 176 0.72 0.60 0.53 0.67 0.62 0.01 2.0 1.99 0.57 -1.10
## se
## Run* 0.24
## conc 0.31
## density 0.04summary(df)## Run conc density
## 10 :16 Min. : 0.04883 Min. :0.0110
## 11 :16 1st Qu.: 0.34180 1st Qu.:0.1978
## 9 :16 Median : 1.17188 Median :0.5265
## 1 :16 Mean : 3.10669 Mean :0.7192
## 4 :16 3rd Qu.: 3.90625 3rd Qu.:1.1705
## 8 :16 Max. :12.50000 Max. :2.0030
## (Other):80# другие статистики можно посчитать примерно так
sapply(df[,-1] , function(x) sd(x))## conc density
## 4.0598653 0.5955726sapply(df[,-1] , function(x) mad(x))## conc density
## 1.5564404 0.6241746sapply(df[,-1] , function(x) mean(sin(x)))## conc density
## 0.2806090 0.5349422# -1 -- это убрать первый столбец, который факторный, а не численныйВизуализации:
#графики
GGally::ggpairs(df,
lower = list(continuous = "points", combo = "box_no_facet", discrete = "facetbar", na =
"na"))## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
ps=ggplot(df,aes(y=density,x=conc,col=Run))+
geom_point(size=3)+
theme_bw()
ps
ps+facet_wrap(.~Run)# из этого графика видно, что зависимость между кол. переменными для каждого уровня Run очень схожая и похожа на sqrt
Тесты:
summary(aov(density ~ conc+Run:conc +Run,df))# дисперсионный анализ подтверждает наличие зависимости, однако от фактора зависимости не выявлено## Df Sum Sq Mean Sq F value Pr(>F)
## conc 1 53.80 53.80 1037.658 <2e-16 ***
## Run 10 0.14 0.01 0.278 0.985
## conc:Run 10 0.15 0.01 0.281 0.985
## Residuals 154 7.98 0.05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cor.test(df$conc,df$density) # тест на корреляцию подтверждает её наличие на весьма высоком уровне##
## Pearson's product-moment correlation
##
## data: df$conc and df$density
## t = 33.635, df = 174, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9081006 0.9482984
## sample estimates:
## cor
## 0.9309673Модели:
# модели
fit1=lm(density~conc,df)
mysummary(fit1)## -----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:
## -----> Ошибка = 0.04701372
##
##
## Call:
## lm(formula = density ~ conc, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30901 -0.19640 -0.03957 0.19498 0.48056
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29488 0.02072 14.23 <2e-16 ***
## conc 0.13657 0.00406 33.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2181 on 174 degrees of freedom
## Multiple R-squared: 0.8667, Adjusted R-squared: 0.8659
## F-statistic: 1131 on 1 and 174 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mdl)
##
## Value p-value Decision
## Global Stat 1.687e+02 0.0000000 Assumptions NOT satisfied!
## Skewness 2.651e+00 0.1034945 Assumptions acceptable.
## Kurtosis 1.235e+01 0.0004405 Assumptions NOT satisfied!
## Link Function 1.537e+02 0.0000000 Assumptions NOT satisfied!
## Heteroscedasticity 5.168e-03 0.9426898 Assumptions acceptable.
##
## -----> БАЗОВЫЕ ГРАФИКИ:
##
## -----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ
##
##
## Shapiro-Wilk normality test
##
## data: mdl$residuals
## W = 0.92034, p-value = 3.249e-08
## [1] 29 45# эта модель статистически значима + R^2 > 80%
# но не прошла тесты на нормальность остатков,
# эксцесс и другие
#
fit2=lm(density~sqrt(conc),df)
mysummary(fit2)## -----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:
## -----> Ошибка = 0.007825747
##
##
## Call:
## lm(formula = density ~ sqrt(conc), data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20009 -0.05060 -0.01110 0.05827 0.30600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.053297 0.011081 -4.81 3.26e-06 ***
## sqrt(conc) 0.550521 0.006287 87.57 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08897 on 174 degrees of freedom
## Multiple R-squared: 0.9778, Adjusted R-squared: 0.9777
## F-statistic: 7668 on 1 and 174 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mdl)
##
## Value p-value Decision
## Global Stat 93.1765 0.0000 Assumptions NOT satisfied!
## Skewness 1.5518 0.2129 Assumptions acceptable.
## Kurtosis 1.1208 0.2897 Assumptions acceptable.
## Link Function 90.3535 0.0000 Assumptions NOT satisfied!
## Heteroscedasticity 0.1504 0.6982 Assumptions acceptable.
##
## -----> БАЗОВЫЕ ГРАФИКИ:
##
## -----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ
##
##
## Shapiro-Wilk normality test
##
## data: mdl$residuals
## W = 0.98126, p-value = 0.01803
## [1] 29 45# эта модель статистически значима + R^2 > 95%
# но не прошла тесты на нормальность остатков,
# надо её ещё улучшить
#
fit3=lm(density~sqrt(conc):Run,df)
mysummary(fit3)## -----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:
## -----> Ошибка = 0.006551772
##
##
## Call:
## lm(formula = density ~ sqrt(conc):Run, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16149 -0.05250 -0.01244 0.06853 0.19823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05330 0.01044 -5.103 9.15e-07 ***
## sqrt(conc):Run10 0.54045 0.01279 42.241 < 2e-16 ***
## sqrt(conc):Run11 0.53912 0.01279 42.136 < 2e-16 ***
## sqrt(conc):Run9 0.53907 0.01279 42.133 < 2e-16 ***
## sqrt(conc):Run1 0.52802 0.01279 41.269 < 2e-16 ***
## sqrt(conc):Run4 0.53422 0.01279 41.754 < 2e-16 ***
## sqrt(conc):Run8 0.54189 0.01279 42.353 < 2e-16 ***
## sqrt(conc):Run5 0.54393 0.01279 42.512 < 2e-16 ***
## sqrt(conc):Run7 0.54876 0.01279 42.890 < 2e-16 ***
## sqrt(conc):Run6 0.55995 0.01279 43.764 < 2e-16 ***
## sqrt(conc):Run2 0.58669 0.01279 45.855 < 2e-16 ***
## sqrt(conc):Run3 0.59363 0.01279 46.397 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08385 on 164 degrees of freedom
## Multiple R-squared: 0.9814, Adjusted R-squared: 0.9802
## F-statistic: 787.7 on 11 and 164 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mdl)
##
## Value p-value Decision
## Global Stat 109.440 0.00000 Assumptions NOT satisfied!
## Skewness 1.589 0.20751 Assumptions acceptable.
## Kurtosis 2.841 0.09187 Assumptions acceptable.
## Link Function 104.200 0.00000 Assumptions NOT satisfied!
## Heteroscedasticity 0.810 0.36813 Assumptions acceptable.
##
## -----> БАЗОВЫЕ ГРАФИКИ:
##
## -----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ
##
##
## Shapiro-Wilk normality test
##
## data: mdl$residuals
## W = 0.9728, p-value = 0.001596
## [1] 45 48# эта модель статистически значима + R^2 > 98%,
# взаимодействие между sqrt(conc) и всеми уровнями Run значимо
# но не прошла тесты на нормальность остатков,
# надо её ещё улучшить
#
fit4=lm(density~sqrt(conc):Run+conc:Run,df)
mysummary(fit4)## -----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:
## -----> Ошибка = 0.002259116
##
##
## Call:
## lm(formula = density ~ sqrt(conc):Run + conc:Run, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09010 -0.03929 -0.01069 0.03669 0.14840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.181979 0.010058 -18.092 < 2e-16 ***
## sqrt(conc):Run10 0.817155 0.026310 31.059 < 2e-16 ***
## sqrt(conc):Run11 0.787265 0.026310 29.923 < 2e-16 ***
## sqrt(conc):Run9 0.787384 0.026310 29.927 < 2e-16 ***
## sqrt(conc):Run1 0.731690 0.026310 27.810 < 2e-16 ***
## sqrt(conc):Run4 0.727955 0.026310 27.668 < 2e-16 ***
## sqrt(conc):Run8 0.771176 0.026310 29.311 < 2e-16 ***
## sqrt(conc):Run5 0.770446 0.026310 29.283 < 2e-16 ***
## sqrt(conc):Run7 0.810580 0.026310 30.809 < 2e-16 ***
## sqrt(conc):Run6 0.816278 0.026310 31.025 < 2e-16 ***
## sqrt(conc):Run2 0.775101 0.026310 29.460 < 2e-16 ***
## sqrt(conc):Run3 0.798293 0.026310 30.342 < 2e-16 ***
## Run10:conc -0.079507 0.008496 -9.358 < 2e-16 ***
## Run11:conc -0.069120 0.008496 -8.136 1.31e-13 ***
## Run9:conc -0.069181 0.008496 -8.143 1.26e-13 ***
## Run1:conc -0.052942 0.008496 -6.231 4.28e-09 ***
## Run4:conc -0.049328 0.008496 -5.806 3.58e-08 ***
## Run8:conc -0.062260 0.008496 -7.328 1.26e-11 ***
## Run5:conc -0.061252 0.008496 -7.210 2.41e-11 ***
## Run7:conc -0.074091 0.008496 -8.721 4.33e-15 ***
## Run6:conc -0.072097 0.008496 -8.486 1.72e-14 ***
## Run2:conc -0.047391 0.008496 -5.578 1.08e-07 ***
## Run3:conc -0.053305 0.008496 -6.274 3.44e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05098 on 153 degrees of freedom
## Multiple R-squared: 0.9936, Adjusted R-squared: 0.9927
## F-statistic: 1079 on 22 and 153 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mdl)
##
## Value p-value Decision
## Global Stat 114.646 0.000000 Assumptions NOT satisfied!
## Skewness 8.517 0.003518 Assumptions NOT satisfied!
## Kurtosis 1.664 0.197009 Assumptions acceptable.
## Link Function 101.690 0.000000 Assumptions NOT satisfied!
## Heteroscedasticity 2.775 0.095774 Assumptions acceptable.
##
## -----> БАЗОВЫЕ ГРАФИКИ:
##
## -----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ
##
##
## Shapiro-Wilk normality test
##
## data: mdl$residuals
## W = 0.95841, p-value = 4.444e-05
## [1] 45 48# модель ещё лучше
fit4=lm(density~sqrt(conc):Run+conc:Run+I(conc^2):Run,df)
mysummary(fit4)## -----> ОБЩАЯ ИНФОРМАЦИЯ О МОДЕЛИ:
## -----> Ошибка = 0.0007179746
##
##
## Call:
## lm(formula = density ~ sqrt(conc):Run + conc:Run + I(conc^2):Run,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.077382 -0.016776 -0.002738 0.012914 0.072618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.081716 0.008416 -9.710 < 2e-16 ***
## sqrt(conc):Run10 0.577542 0.030742 18.786 < 2e-16 ***
## sqrt(conc):Run11 0.530702 0.030742 17.263 < 2e-16 ***
## sqrt(conc):Run9 0.501568 0.030742 16.315 < 2e-16 ***
## sqrt(conc):Run1 0.448894 0.030742 14.602 < 2e-16 ***
## sqrt(conc):Run4 0.406788 0.030742 13.232 < 2e-16 ***
## sqrt(conc):Run8 0.469522 0.030742 15.273 < 2e-16 ***
## sqrt(conc):Run5 0.462764 0.030742 15.053 < 2e-16 ***
## sqrt(conc):Run7 0.584415 0.030742 19.010 < 2e-16 ***
## sqrt(conc):Run6 0.569106 0.030742 18.512 < 2e-16 ***
## sqrt(conc):Run2 0.436930 0.030742 14.213 < 2e-16 ***
## sqrt(conc):Run3 0.511337 0.030742 16.633 < 2e-16 ***
## Run10:conc 0.035732 0.018561 1.925 0.05621 .
## Run11:conc 0.058940 0.018561 3.176 0.00183 **
## Run9:conc 0.081009 0.018561 4.365 2.44e-05 ***
## Run1:conc 0.094964 0.018561 5.116 9.93e-07 ***
## Run4:conc 0.127604 0.018561 6.875 1.81e-10 ***
## Run8:conc 0.099911 0.018561 5.383 2.96e-07 ***
## Run5:conc 0.105479 0.018561 5.683 7.24e-08 ***
## Run7:conc 0.030974 0.018561 1.669 0.09736 .
## Run6:conc 0.048859 0.018561 2.632 0.00942 **
## Run2:conc 0.142404 0.018561 7.672 2.46e-12 ***
## Run3:conc 0.097747 0.018561 5.266 5.04e-07 ***
## Run10:I(conc^2) -0.004537 0.000903 -5.024 1.50e-06 ***
## Run11:I(conc^2) -0.005214 0.000903 -5.774 4.68e-08 ***
## Run9:I(conc^2) -0.006384 0.000903 -7.070 6.46e-11 ***
## Run1:I(conc^2) -0.006263 0.000903 -6.936 1.31e-10 ***
## Run4:I(conc^2) -0.007798 0.000903 -8.635 1.09e-14 ***
## Run8:I(conc^2) -0.007018 0.000903 -7.771 1.43e-12 ***
## Run5:I(conc^2) -0.007259 0.000903 -8.038 3.22e-13 ***
## Run7:I(conc^2) -0.003999 0.000903 -4.428 1.88e-05 ***
## Run6:I(conc^2) -0.004839 0.000903 -5.358 3.31e-07 ***
## Run2:I(conc^2) -0.008478 0.000903 -9.388 < 2e-16 ***
## Run3:I(conc^2) -0.006430 0.000903 -7.120 4.93e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02983 on 142 degrees of freedom
## Multiple R-squared: 0.998, Adjusted R-squared: 0.9975
## F-statistic: 2109 on 33 and 142 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma::gvlma(x = mdl)
##
## Value p-value Decision
## Global Stat 79.5500 2.220e-16 Assumptions NOT satisfied!
## Skewness 4.0543 4.406e-02 Assumptions NOT satisfied!
## Kurtosis 1.8792 1.704e-01 Assumptions acceptable.
## Link Function 73.0629 0.000e+00 Assumptions NOT satisfied!
## Heteroscedasticity 0.5536 4.568e-01 Assumptions acceptable.
##
## -----> БАЗОВЫЕ ГРАФИКИ:
##
## -----> ТЕСТ НА НОРМАЛЬНОСТЬ РАСПРЕДЕЛЕНИЯ ОСТАТКОВ
##
##
## Shapiro-Wilk normality test
##
## data: mdl$residuals
## W = 0.97446, p-value = 0.002523
## [1] 45 46# модель ещё лучше, лучше этой линейную модель уже вряд ли сделаешь
# в документации есть ещё такие нелинейные модели, но у них ошибка больше даже чем у fit3
fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
data = DNase, subset = Run == 1)
## compare with a four-parameter logistic
fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
data = DNase, subset = Run == 1)
summary(fm2)##
## Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A -0.007897 0.017200 -0.459 0.654
## B 2.377239 0.109516 21.707 5.35e-11 ***
## xmid 1.507403 0.102080 14.767 4.65e-09 ***
## scal 1.062579 0.056996 18.643 3.16e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01981 on 12 degrees of freedom
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 2.517e-07anova(fm1, fm2)## Analysis of Variance Table
##
## Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
## Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
## 1 13 0.0047896
## 2 12 0.0047073 1 8.2314e-05 0.2098 0.6551(predict(fm1,df)-df$density)^2 %>% mean()## [1] 0.003935952(predict(fm2,df)-df$density)^2 %>% mean()## [1] 0.003984377