devbox
我家自动化
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

人工晶状体计算——人工智能算法（R语言）_用r语言做人工智能

作者：我家自动化 | 2024-03-22 19:13:33

赞

踩

用r语言做人工智能

人工晶状体计算——人工智能算法（R语言）

1. 准备数据

2. 建立模型

2.1 方法1

2.2 方法2

1. 准备数据

准备数据Data.xlsx，示例如图

 

AgeALACDK1K2WTWLTRefIOLPower
68.0022.682.2142.4442.7511.204.82-1.1326.50
62.0023.793.4243.9345.5111.704.49-0.1119.50
62.0023.823.1244.0844.5312.004.90-0.8621.00
68.0022.562.3042.6543.1311.004.65-0.3725.50
73.0023.502.3843.4343.8411.404.85-0.2621.50
74.0023.092.9143.8644.7411.904.45-0.5323.00
73.0024.613.5741.7042.2311.903.49-0.8021.00
64.0022.472.8946.6348.0111.104.41-0.5221.00
63.0021.352.6743.8944.5312.204.68-0.5829.50

2. 建立模型

2.1 方法1

直接建立Age,AL,ACD,K1,K2,WTW,LT,Ref与IOLPower之间的关系

将使用多种人工智能模型，最后将表现最优的三种叠加在一起。

开始先看一下各自变量之间的相关性：

可以看到自变量之间还是有一定相关性的，比如LT(Lens Thickness)与Age正相关，ACD与AL正相关等。

  1. rm(list=ls())
  2. # load libraries
  3. library(caret)
  4. library(corrplot)
  5. library(tidyverse)
  6. library(xlsx)
  7. # Split out validation dataset
  8. set.seed(123)
  9. df <- read.xlsx("./Data.xlsx",sheetName ="Sheet1",stringsAsFactors = FALSE)
  10. str(df)
  11. names(df)<-c("Age","AL","ACD","K1","K2","WTW","LT","IOLPower","Ref")
  12. cor(df[c("Age","AL","ACD","K1","K2","WTW","LT")])
  13. pairs(df[c("Age","AL","ACD","K1","K2","WTW","LT")])
  14. # more informative scatterplot matrix
  15. library(psych)
  16. pairs.panels(df[c("Age","AL","ACD","K1","K2","WTW","LT")])
  17.  
  18. train.idx <- sample(1:nrow(df), 0.7*nrow(df))
  19. train.df <- as_tibble(df[train.idx,])
  20. test.df <- as_tibble(df[-train.idx,])
  21. train_x = train.df %>% select(-'IOLPower')
  22. train_y = train.df[,'IOLPower']
  23. test_x = test.df %>% select(-'IOLPower')
  24. test_y = test.df[,'IOLPower']
  25. # Run algorithms using 10-fold cross validation
  26. control <- trainControl(method="repeatedcv", number=10, repeats=3, search="random")
  27. metric <- "RMSE"
  28. # lm
  29. set.seed(123)
  30. fit_lm <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data=train.df, method="lm",
  31. metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  32. fit_lm
  33. # SVM
  34. set.seed(123)
  35. fit_svm <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df,
  36. method="svmRadial", metric=metric, preProc=c("center", "scale"), tuneLength=15,
  37. trControl=control)
  38. fit_svm
  39.  
  40. # ANNs
  41. fit_nnet <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df, method="nnet",
  42.                   metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  43. fit_nnet
  44. if(FALSE)
  45. {
  46.   pred_Power <-predict(fit_nnet,test_x)
  47.   print(pred_Power)
  48.   err<- test_y$IOLPower-pred_Power
  49.   print(err)
  50.   plot(err)
  51. }
  52.  
  53. # random forest
  54. set.seed(123)
  55. fit_rf <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.df, method="rf",
  56. metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  57. fit_rf
  58. # CART
  59. set.seed(123)
  60. fit_cart <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df, method="rpart",
  61. metric=metric,preProc=c("center", "scale"),tuneLength=15, trControl=control)
  62. fit_cart
  63. # kNN
  64. set.seed(123)
  65. fit_knn <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df, method="knn",
  66. metric=metric, preProc=c("center", "scale"),tuneLength=15, trControl=control)
  67. fit_knn
  68. #MARS
  69. set.seed(123)
  70. fit_MARS<- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df,
  71. method="earth", metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  72. fit_MARS
  73. # GBM
  74. set.seed(123)
  75. fit_gbm <- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.df, method =
  76. "gbm",metric=metric,preProc=c("center", "scale"), tuneLength=15, trControl=control)
  77. fit_gbm
  78. #cubist
  79. set.seed(123)
  80. trControl=control
  81. fit_cubist<- train(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.df, method =
  82.                      "cubist",metric=metric,preProc=c("center", "scale"), tuneLength=15, trControl=control)
  83. fit_cubist
  84. # Compare algorithms
  85. results <- resamples(list(LM = fit_lm, CART = fit_cart, KNN= fit_knn, RF = fit_rf, SVR = fit_svm,
  86. GBM = fit_gbm, MARS= fit_MARS, Cubist = fit_cubist,ANNs=fit_nnet))
  87. summary(results)
  88. dotplot(results)
  89. # Stacking
  90. library(caretEnsemble)
  91. set.seed(123)
  92.  
  93. allModels <- caretList(IOLPower ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.df,trControl
  94. = trainControl("repeatedcv", number = 10, repeats = 3, search="random",verboseIter = T),metric =
  95. "RMSE", tuneList = list(cubist = caretModelSpec(method = "cubist", tuneLength = 15, preProc =
  96. c("center","scale")),svmRadial = caretModelSpec(method = "svmRadial", tuneLength = 15, preProc
  97. = c("center","scale")),Earth = caretModelSpec(method = "earth", tuneLength = 15, preProc =
  98. c("center","scale"))))
  99.  
  100. bwplot(resamples(allModels), metric = "RMSE")
  101. summary(allModels)
  102. # Stack with lm
  103. set.seed(123)
  104. stack <- caretStack(allModels, metric = "RMSE", method = "lm", tuneLength = 10,trControl =
  105. trainControl("repeatedcv", number = 5, repeats = 3, verboseIter = T))
  106. summary(stack)
  107. print(stack)
  108.  
  109. if(TRUE){
  110.   ggplot(varImp(fit_lm))
  111.   pred_Power <-predict(stack,test_x)
  112.   print(pred_Power)
  113.   plot(pred_Power)
  114.   pred_err <- test_y$IOLPower - pred_Power
  115.   print(pred_err)
  116.   summary(pred_err)
  117.   plot(pred_err)
  118. }

运行结果误差还是很小的。

2.2 方法2

结合光学公式与人工智能（此法源自中山的一篇论文）

光学人工晶体计算公式如下

Power =1336 / (AL- ELP) - (1336/(1336 /(1000 /(1000 /Ref-12)+K)-ELP))

临床实践中最难的就是预测ELP(Effective Lens Position，有效人工晶状体位置，个人觉得"等效人工晶状体位置"这个翻译可能更好)，以下代码中Ratio指

Ratio =(ELP-ACD)/LT

下面代码的各个模型都是为了预测这个Ratio，最后求得ELP，然后代入上述光学计算公式中得到最后的人工晶状体度数。

  1. if (TRUE){
  2.   rm(list=ls())
  3.   library(xlsx)
  4.   data <- read.xlsx("./Data.xlsx",sheetName ="Sheet1",stringsAsFactors = FALSE)
  5.   f<-function(ELP,AL,K,Ref,P)
  6.   {1336 / (AL- ELP) - (1336/(1336 /(1000 /(1000 /Ref-12)+K)-ELP))-P}
  7.   root<-c()
  8.   
  9.   for(i in 1:nrow(data)){
  10.     root[i]<-uniroot(f,c(0,20),AL=data$AL[i],K=(data$K1[i]+data$K2[i])/2,R=data$Ref[i],P=data$IOLPower[i])$root
  11.   }
  12.   root
  13.   data$ELP <-root
  14.   Ratio <- (data$ELP-data$ACD)/data$LT
  15.   Ratio
  16.   data$Ratio<-Ratio
  17.   PowerCalc<-function(ELP,AL,K,Ref)
  18.   {1336 / (AL- ELP) - (1336/(1336 /(1000 /(1000 /Ref-12)+K)-ELP))}
  19.   P<-PowerCalc(data$ELP,data$AL,(data$K1+data$K2)/2,data$Ref)
  20.   P
  21.   options(scipen = 100) #小数点后100位不使用科学计数法
  22.   err <- data$IOLPower - P
  23.   summary(err)
  24.   plot(err)
  25. }
  26.  
  27.  
  28. library(caret)
  29. library(corrplot)
  30. library(tidyverse)
  31. library(xlsx)
  32. # Split out validation dataset
  33. set.seed(123)
  34. cor(data[c("Age","AL","ACD","K1","K2","WTW","LT")])
  35. pairs(data[c("Age","AL","ACD","K1","K2","WTW","LT")])
  36. # more informative scatterplot matrix
  37. library(psych)
  38. pairs.panels(data[c("Age","AL","ACD","K1","K2","WTW","LT")])
  39.  
  40. train.idx <- sample(1:nrow(data), 0.7*nrow(data))
  41. train.data <- as_tibble(data[train.idx,])
  42. test.data <- as_tibble(data[-train.idx,])
  43. train_x = train.data %>% select(-'Ratio')
  44. train_y = train.data[,'Ratio']
  45. test_x = test.data %>% select(-'Ratio')
  46. test_y = test.data[,'Ratio']
  47. # Run algorithms using 10-fold cross validation
  48. control <- trainControl(method="repeatedcv", number=10, repeats=3, search="random")
  49. metric <- "RMSE"
  50. # lm
  51. set.seed(123)
  52. fit_lm <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data=train.data, method="lm",
  53.                 metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  54. fit_lm
  55. # SVM
  56. set.seed(123)
  57. fit_svm <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data,
  58.                  method="svmRadial", metric=metric, preProc=c("center", "scale"), tuneLength=15,
  59.                  trControl=control)
  60. fit_svm
  61. plot(fit_svm)
  62. # ANNs
  63. fit_nnet <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data, method="nnet",
  64.                   metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  65. fit_nnet
  66. if(TRUE)
  67. {
  68.   pred <-predict(fit_nnet,test_x)
  69.   print(pred)
  70.   err<- test_y$Ratio-pred
  71.   print(err)
  72.   plot(err)
  73. }
  74.  
  75. # random forest
  76. set.seed(123)
  77. fit_rf <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.data, method="rf",
  78.                 metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  79. fit_rf
  80. # CART
  81. set.seed(123)
  82. fit_cart <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data, method="rpart",
  83.                   metric=metric,preProc=c("center", "scale"),tuneLength=15, trControl=control)
  84. fit_cart
  85. # kNN
  86. set.seed(123)
  87. fit_knn <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data, method="knn",
  88.                  metric=metric, preProc=c("center", "scale"),tuneLength=15, trControl=control)
  89. fit_knn
  90. #MARS
  91. set.seed(123)
  92. fit_MARS<- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data,
  93.                  method="earth", metric=metric, preProc=c("center", "scale"), tuneLength=15, trControl=control)
  94. fit_MARS
  95. # GBM
  96. set.seed(123)
  97. fit_gbm <- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.data, method =
  98.                    "gbm",metric=metric,preProc=c("center", "scale"), tuneLength=15, trControl=control)
  99. fit_gbm
  100. #cubist
  101. set.seed(123)
  102. trControl=control
  103. fit_cubist<- train(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref, data = train.data, method =
  104.                      "cubist",metric=metric,preProc=c("center", "scale"), tuneLength=15, trControl=control)
  105. fit_cubist
  106. # Compare algorithms
  107. results <- resamples(list(LM = fit_lm, CART = fit_cart, KNN= fit_knn, RF = fit_rf, SVR = fit_svm,
  108.                           GBM = fit_gbm, MARS= fit_MARS, Cubist = fit_cubist,ANNs=fit_nnet))
  109. summary(results)
  110. dotplot(results)
  111. # Stacking
  112. library(caretEnsemble)
  113. set.seed(123)
  114.  
  115. allModels <- caretList(Ratio ~ Age+AL+ACD+K1+K2+WTW+LT+Ref,data=train.data,trControl
  116.                        = trainControl("repeatedcv", number = 10, repeats = 3, search="random",verboseIter = T),metric =
  117.                          "RMSE", tuneList = list(cubist = caretModelSpec(method = "cubist", tuneLength = 15, preProc =
  118.                                                                            c("center","scale")),svmRadial = caretModelSpec(method = "svmRadial", tuneLength = 15, preProc
  119.                                                                                                                            = c("center","scale")),Earth = caretModelSpec(method = "earth", tuneLength = 15, preProc =
  120.                                                                                                                                                                            c("center","scale"))))
  121.  
  122. bwplot(resamples(allModels), metric = "RMSE")
  123. summary(allModels)
  124. # Stack with lm
  125. set.seed(123)
  126. stack <- caretStack(allModels, metric = "RMSE", method = "lm", tuneLength = 10,trControl =
  127.                       trainControl("repeatedcv", number = 5, repeats = 3, verboseIter = T))
  128. summary(stack)
  129. print(stack)
  130.  
  131. if(TRUE){
  132.   ggplot(varImp(fit_nnet))
  133.   pred_Ratio <-predict(stack,test_x)
  134.   print(pred_Ratio)
  135.   plot(pred_Ratio)
  136.   pred_err1 <- test_y$Ratio - pred_Ratio
  137.   print(pred_err1)
  138.   summary(pred_err1)
  139.   plot(pred_err1)
  140.   ELPCalc<-function(Ratio,LT,ACD){Ratio*LT+ACD}
  141.   ELP = ELPCalc(pred_Ratio,test_x$LT,test_x$ACD)
  142.   Pred_Power<-PowerCalc(ELP,test_x$AL,(test_x$K1+test_x$K2)/2,test_x$Ref)
  143.   Pred_err2<-test_x$IOLPower - Pred_Power
  144.   summary(Pred_err2)
  145.   plot(Pred_err2)
  146. }

复现论文也算勉强成功，但是样本量太少，暂无法比较两种方法的优劣。

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/我家自动化/article/detail/290335
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号