运行一个 for 循环，其中模拟数据然后平均

问题描述

我想运行以下代码 100 次：

require(phytools)
require(RRphylo)

linbm <- rowMeans(fastBM(lattree,a = spreadlind[1,1],bounds=c(min(spreadlind),max(spreadlind)),nsim = 100))
RRlinbm <- RRphylo(lattree,linbm)
trendlinbm <- search.trend(RRlinbm,linbm,foldername=tempdir())

spreadlindbm <- trendlinbm$pbt[,1:2]
spreadlindbm1 <- spreadlindbm[order(spreadlindbm$age),]
spreadlindbm2 <- (max(spreadlindbm1$phenotype)-min(spreadlindbm1$phenotype))
spreadlindbm3 <- spreadlindbm1[which(spreadlind1$age < ((max(spreadlindbm1$age)-min(spreadlindbm1$age))/2+min(spreadlindbm1))),]
spreadlindbm4 <- (max(spreadlindbm3$phenotype)-min(spreadlindbm3$phenotype))
spreadlindbm5 <- spreadlindbm2/spreadlindbm4

spreads <- spreadlind5/spreadlindbm5

这很好用（尽管有一些警告）并且可以完成我想要它做的事情。但是，当我这样做时：

linbm <- list()
RRlinbm <- list()
trendlinbm <- list()
spreadlindbm <- list()
spreadlindbm1 <- list()
spreadlindbm2 <- list()
spreadlindbm3 <- list()
spreadlindbm4 <- list()
spreadlindbm5 <- list()
spreads <- list()

for (i in 1:100) {
linbm[[i]] <- rowMeans(fastBM(lattree,1][[i]],bounds=c(min(spreadlind[[i]]),max(spreadlind[[i]])),nsim = 100))
RRlinbm[[i]] <- RRphylo(lattree,linbm[[i]])
trendlinbm[[i]] <- search.trend(RRlinbm[[i]],linbm[[i]],foldername=tempdir())

spreadlindbm[[i]] <- trendlinbm$pbt[,1:2][[i]]
spreadlindbm1[[i]] <- spreadlindbm[[i]][order(spreadlindbm$age[[i]]),]
spreadlindbm2[[i]] <- (max(spreadlindbm1$phenotype[[i]])-min(spreadlindbm1$phenotype[[i]]))
spreadlindbm3[[i]] <- spreadlindbm1[[i]][which(spreadlind1$age[[i]] < ((max(spreadlindbm1$age[[i]])-min(spreadlindbm1$age[[i]]))/2+min(spreadlindbm1[[i]]))),]
spreadlindbm4[[i]] <- (max(spreadlindbm3$phenotype[[i]])-min(spreadlindbm3$phenotype[[i]]))
spreadlindbm5[[i]] <- spreadlindbm2[[i]]/spreadlindbm4[[i]]

spreads[[i]] <- spreadlind5[[i]]/spreadlindbm5[[i]]
}

我收到此错误：

spreadlindbm[[i]] 错误：下标越界

似乎只有 1 次模拟 linbm 而不是 100，这就是我收到错误的原因。

我认为这与我似乎无法执行 lattree[[i]] 的事实有关，因为如果我这样做，出于某种原因，它会调用一个数据集，显示哪些节点连接到哪些节点，而不是树本身，并且因此我最终得到一个错误，指出该对象不再属于 phylo 类。

因此我不确定如何解决这个问题。

循环中的警告与原始脚本相同。

任何帮助将不胜感激！

谢谢，

卡罗莱纳州

这是数据：

> dput(spreadlind)
structure(list(y.multi = c(0.00155347055674816,0.161613111788121,0.782198041620036,1.53091468559125,1.5374150126118,1.95106361585544,2.84877753353601,3.30856549415935,3.66347992385696,3.61603383156884,4.14173027008837,4.81866865201557,5.96045394613843,3.62927975265028,3.63693964926282,6.36794280876027,6.31025739368824,3.68593745526834,3.97891742720543,6.93648438239245,0.494718563277189,0.582320520627804,1.80998664230687,2.07126271026748,2.09882261692225,2.04819403524499,2.67566519385064,2.72232617371275,0.51140775510139,3.54745844047975,1.25074485672904,1.55010916708204,1.2871994431796,1.32085767556071,1.37033008644747,0.515214782210954,0.885611593749094,0.927624990707522,1.2353386216543,1.20393998975362,1.05078386017598,1.53796097273876,1.43610602382877,0.571743119736813,3.59952827678168,0.571885350269405,1.2853590224863,0.579301896064088,2.86926984203778,0.587708152457224,1.82659625228065,3.60766146697518,0.494882679677509,0.506429821730035,0.522518109651434,1.90891540289338,2.3451011743209,2.33727446999219,0.523615252037944,2.36609248926363,2.59202193074751,2.66474569109868,1.72333855229808,2.4163690782667,1.73287568525806,2.49842570272719,1.80664555049829,2.38385429265236,2.73704621987446,0.506784376709435,0.507707989713869,2.2820288381412,1.92420878906418,1.94303773985105,2.55650284133446,2.23982669071852,2.62370104646633,2.74849048160226,2.7516366913285,0.540677960556539,0.565447800474045,2.97105772656143,1.26383072332201,0.543457139827625,3.58939916599798,1.11500285637845,1.14441613040443,1.19209724885874,0.783385958851092,1.91325256635718,3.04797028869207,1.92396908318128,4.08905782084652,1.97222419679872,3.67244745727709,2.07774303591261,3.90466176531347,2.91480361864387,0.798118906868094,1.45926213747926,3.52510770252059,1.58839858730842,1.79707280379355,2.47505145152866,2.42499608715655,1.59253486654689,2.590738827913,1.64370161579295,3.81646886358259,3.82520657528986,1.33563102388181,1.43009331011062,2.56134735686846,1.49433262087347,3.27754161801841,1.50755033198619,2.45008826727309,2.43835482679194,1.33914528685374,2.67199581109489,1.68819329149765,2.06661520323385,1.69782065598436,2.37622959857406,2.3078393520338,2.64428969434044,2.31132541432424,2.79957533334395,2.68374060836952),age = c(1.0001,22.8060363126218,54.5976175658671,69.753316771761,69.7533167717499,80.414170921278,88.2586282621069,93.7004406841295,93.7004406841094,90.9184836136638,95.0615012270222,101.543733610294,101.543733610292,91.4587431450015,92.8185903007135,99.6528647376402,99.6528647376378,92.5908102945527,101.011824568968,101.011824568975,18.259071508084,25.993913616012,59.8238331049015,68.743553809163,68.7435538091594,86.2934244425716,102.57347329813,102.573473298177,21.1774077435026,96.7763651286279,80.3066760159123,98.0397518161518,86.9698144837543,98.0397518161568,98.0397518161501,22.741886784047,42.9286275525969,52.0750257961256,66.3690315373348,66.3690315373326,56.6487478575139,75.5154297808487,75.5154297807598,28.0487885610707,95.7021705494929,28.3608373184629,95.3901217920203,30.6460556878806,95.3901217919969,33.5941512838692,92.4420261960504,92.4420261960636,18.6204905115437,20.3679013654434,22.560804157767,88.2577872001976,101.080471597327,101.080471597338,23.1614011949379,84.3358417514739,100.479874560188,100.479874560183,43.8197893493495,79.8214864058532,45.7871594004954,77.8541163546938,51.6615048191536,71.979770936008,71.9797709360253,21.0176870990147,22.0184256013505,99.4647314842312,67.2132181179478,71.8437889797186,94.8341606223199,94.83416062225,97.4250876440769,99.4647314840716,99.4647314840717,29.1448747815453,35.5461381027468,94.05802155308,94.0580215530736,31.2193215751653,100.459284874452,84.0675074329436,100.459284874353,100.45928487436,27.2065400456165,44.9626265916094,92.3960105252684,46.5900338274134,90.8408601707908,50.2586889363923,87.2788979028246,56.4058856487657,81.1317011904194,81.1317011904359,30.8212285638619,48.4056327901052,88.8880148481789,55.1125261701695,62.7514247521627,81.0334816994758,81.033481699512,55.9884815439191,88.8293630359064,58.930241329349,86.066774970786,86.0667749707636,49.3396177771712,54.4444816307805,70.3696256347752,60.2061611560702,70.2835027623113,61.0861700776571,70.4929740447419,70.4208129948606,50.5926362075624,69.1166072043804,62.7279938432217,69.116607204367,63.6382587409568,68.2063423066373,67.4018609248549,68.1881923860238,67.4823724506412,68.1491655977022,68.1463399182352
)),class = "data.frame",row.names = c("67","68","69","Meleagris_gallopavo","Lagopus_lagopus","70","71","Branta_canadensis","Branta_bernicla","72","73","Anas_crecca","Anas_clypeata","74","75","Mergus_merganser","Bucephala_islandica","76","Aix_sponsa","Tadorna_tadorna","77","78","79","Caprimulgus_climacurus","Caprimulgus_europaeus","80","Apus_apus","Tachymarptis_melba","81","82","Podiceps_cristatus","83","84","Streptopelia_decaocto","85","Columba_livia","Columba_palumbus","86","87","88","Fulica_atra","gallinula_chloropus","89","Rallus_aquaticus","Crex_crex","90","Gavia_stellata","91","Fulmarus_glacialis","92","Platalea_leucorodia","93","Fregata_aquila","Phalacrocorax_carbo","94","95","96","97","Falco_peregrinus","Falco_columbarius","98","99","Hirundo_rustica","Delichon_urbicum","100","Amazona_aestiva","101","Psittacula_krameri","102","Trichoglossus_haematodus","Melopsittacus_undulatus","103","104","Sarcoramphus_papa","105","106","Accipiter_nisus","Accipiter_gentilis","107","Aquila_pomarina","Aquila_clanga","108","109","Anthracoceros_coronatus","Upupa_epops","110","Alcedo_atthis","111","Picus_viridis","Dendrocopos_major","112","113","Pluvialis_apricaria","114","Haematopus_ostralegus","115","Vanellus_vanellus","116","Charadrius_dubius","Charadrius_hiaticula","117","118","Numenius_arquata","119","120","Philomachus_pugnax","Calidris_alpina","121","Actitis_hypoleucos","122","Scolopax_rusticola","gallinago_gallinago","123","124","Fratercula_arctica","125","Uria_aalge","126","Alca_torda","Alle_alle","127","Chlidonias_leucopterus","128","Rissa_tridactyla","129","Larus_ridibundus","130","Larus_glaucoides","131","Larus_marinus","Larus_argentatus"
))

还有我的树：

> dput(lattree)
structure(list(edge = structure(c(67L,68L,69L,70L,71L,72L,73L,74L,75L,76L,67L,77L,78L,79L,80L,81L,82L,83L,84L,85L,86L,87L,88L,89L,90L,91L,92L,93L,94L,95L,96L,97L,98L,99L,100L,101L,102L,103L,104L,105L,106L,107L,108L,109L,110L,111L,112L,113L,114L,115L,116L,117L,118L,119L,120L,121L,122L,123L,124L,125L,126L,127L,128L,129L,1L,2L,3L,4L,9L,10L,5L,6L,7L,8L,63L,64L,65L,66L,20L,21L,22L,23L,11L,12L,13L,14L,19L,15L,18L,16L,17L,61L,62L,59L,60L,58L,57L,55L,56L,45L,46L,47L,48L,49L,50L,51L,52L,53L,54L,44L,43L,42L,40L,41L,34L,35L,36L,39L,37L,38L,30L,33L,31L,32L,24L,29L,25L,28L,26L,27L),.Dim = c(128L,2L)),edge.length = c(21.8060363126218,31.7915812532453,15.155699205894,15.1556992058829,57.6081346086562,7.84445734082895,5.44181242202259,5.44181242200252,10.5043126923859,4.14301761335842,6.48223238327186,6.48223238327021,0.540259531337661,1.35984715571197,6.8342744369267,6.83427443692435,1.13206714955123,8.42101427441535,8.42101427442253,17.259071508084,7.73484210792803,33.8299194888895,8.91972070426153,8.91972070425791,60.2995108265596,16.2800488555586,16.2800488556051,2.91833623541861,75.5989573851252,59.1292682724097,17.7330758002395,6.66313846784204,11.0699373324024,11.0699373323958,1.56447904054438,20.1867407685499,9.14639824352872,14.2940057412092,14.294005741207,13.720120304917,18.8666819233348,18.8666819232459,5.30690177702373,67.6533819884221,0.312048757392168,67.0292844735574,2.28521836941774,64.7440661041163,2.94809559598857,58.8478749121812,58.8478749121944,0.361419003459647,1.74741085389973,2.19290279232366,65.6969830424305,12.8226843971294,12.8226843971402,0.600597037170883,61.174440556536,16.1440328087146,16.1440328087093,20.6583881544116,36.0016970565037,1.96737005114587,32.0669569541984,5.87434541865818,20.3182661168544,20.3182661168717,0.649785733571314,1.00073850233577,77.4463058828807,45.1947925165973,4.63057086177076,22.9903716426013,22.9903716425314,30.2118695261291,2.03964383999477,2.03964383999479,8.12718768253063,6.40126332120143,58.5118834503332,58.5118834503268,2.07444679362001,69.2399632992864,52.8481858577782,16.3917774414099,16.3917774414166,8.58604953407289,17.7560865459928,47.4333839336591,1.627407235804,44.2508263433775,3.6686551089789,37.0202089664324,6.14719671237347,24.7258155416537,24.7258155416702,3.61468851824533,17.5844042262434,40.4823820580737,6.70689338006431,7.6388985819932,18.2820569473131,18.2820569473492,0.875955373749606,32.8408814919872,2.94175978542986,27.136533641437,27.1365336414146,18.5183892133093,5.10486385360935,15.9251440039946,5.76167952528963,10.0773416062411,0.880008921586962,9.40680396708479,9.33464291720349,1.25301843039117,18.523970996818,12.1353576356593,6.38861336114526,0.910264897735043,4.56808356568051,3.76360218389816,0.786331461168885,0.0805115257862825,0.666793147061034,0.663967467593978),Nnode = 63L,root.edge = 0,tip.label = c("Meleagris_gallopavo","Larus_argentatus","Tachymarptis_melba")),class = "phylo",order = "cladewise")

解决方法

所以经过多次尝试，这个解决方案似乎有效：

linbm <- list()
linbm2 <- list()
RRlinbm <- list()
trendlinbm <- list()
spreadlindbm <- list()
spreadlindbm1 <- list()
spreadlindbm2 <- list()
spreadlindbm3 <- list()
spreadlindbm4 <- list()
spreadlindbm5 <- list()
spreads <- list()
lintrees <- rep.multiPhylo(lattree,100)
alin <- rep(spreadlind[1,1],100)
minlin <- rep(min(spreadlind),100)
maxlin <- rep(max(spreadlind),100)

for (i in 1:100) {
  linbm[[i]] <- fastBM(lintrees[[i]],a = alin[[i]],bounds = c(minlin[[i]],maxlin[[i]]),nsim = 100)
}

for (i in 1:100){  
  linbm2[[i]] <- rowMeans(linbm[[i]])
}  

for (i in 1:100){
  RRlinbm[[i]] <- RRphylo(lintrees[[i]],linbm2[[i]])
}  

for (i in 1:100){
  trendlinbm[[i]] <- search.trend(RRlinbm[[i]],linbm2[[i]],foldername=tempdir())
}

for (i in 1:100){  
  spreadlindbm[[i]] <- trendlinbm[[i]]$pbt
}  

for (i in 1:100){    
  spreadlindbm1[[i]] <- spreadlindbm[[i]][order(spreadlindbm[[i]]$age),]
}

for (i in 1:100){  
  spreadlindbm2[[i]] <- (max(spreadlindbm1[[i]]$phenotype)-min(spreadlindbm1[[i]]$phenotype))
}  

for (i in 1:100){    
  spreadlindbm3[[i]] <- spreadlindbm1[[i]][which(spreadlind1$age < ((max(spreadlindbm1[[i]]$age)-min(spreadlindbm1[[i]]$age))/2+min(spreadlindbm1[[i]]))),]
}

for (i in 1:100){  
  spreadlindbm4[[i]] <- (max(spreadlindbm3[[i]]$phenotype)-min(spreadlindbm3[[i]]$phenotype))
}

for (i in 1:100){  
  spreadlindbm5[[i]] <- spreadlindbm2[[i]]/spreadlindbm4[[i]]
}

for (i in 1:100){    
  spreads[[i]] <- spreadlind5/spreadlindbm5[[i]]
}

虽然它绝对不优雅，所以如果有人有更好的解决方案，请发布它，我会接受它作为解决方案:)