时间戳

首页 > 解决方案 > 用循环计算摄入能量

r - 用循环计算摄入能量

问题描述

我正在尝试运行一个旧的同事脚本，我希望有人可以帮助告诉我他在这段代码段中到底做了什么。在脚本的前面，我们计算了几种猎物的摄入率，现在看来我们正在根据独特的位置对它们进行分组。此后的代码部分要求有 41 行（完整数据集中每个唯一位置 1 行）。我相信代码会根据纬度对数据进行子集化，然后添加一个“alpha”列。我遇到的主要问题是这条线在计算什么：x= x + d$Intakerate_kjday[j]*d$alpha[j]。对于有多个猎物的位置（profit.fall.38.4.959），这段代码是否将“intakerate_kjday”和“alpha”相加，然后将它们相乘？执行代码时，我收到错误 Error in

 `[<-.data.frame`(`*tmp*`, k, , value = c("2", "Bishop's Head",  : replacement has 6 items, need 7

对于他试图计算的内容和潜在的解决方法，我真的很感激。谢谢你。

dput(profit)
structure(list(Sample.ID = structure(c(5L, 19L, 27L, 28L, 30L, 
38L, 12L, 62L, 49L, 29L, 25L, 17L, 61L, 67L, 27L, 26L, 32L, 9L, 
47L, 45L, 5L, 26L, 27L, 44L, 45L, 4L, 1L, 43L, 19L, 35L), .Label = c("Barren Island Mud 1", 
"BH High 1", "BH High 2", "BH Low 1", "BH Low 2", "BH Low 3", 
"BHH 1 C", "BHH 2 E", "BHL 1 E", "BHL 2", "BHL 3 (B)", "BHM 1 C", 
"BI High 1", "BI Low 1", "BI Low 2C", "BI Low 3", "BI Mud", "BIHI High B", 
"BIL1 (low) E", "BIL1E", "BIL2 E", "BIL2E", "BW Fresh 1", "BW Fresh 2", 
"BW High 1", "BW High 2", "BW High 5", "BW Low 3", "BW Money Stump", 
"BW Mud 1", "BW SAV 1", "BW SAV 2", "BWH 1 D", "BWH 2", "BWH 3", 
"BWH 5", "BWL 1", "BWL 2", "BWL 3", "BWM 1", "BWMS D", "EN High 2", 
"EN High 4", "EN High 5", "EN Low 1", "EN Low 2", "EN Mud 2", 
"ENH3 A High", "ENH4 A High", "ENH5 A High", "ENL1 Low E", "ENM1 A Mud", 
"ENS1 SAV", "ENS2 SAV 2C", "ENS3 SAV 3E", "High 3C", "MWP 29 Low 1", 
"MWP 30 Mud 1", "MWP 31 Low 2", "MWP 32 Mud 2", "MWP 33 Low 3", 
"MWP 34 Low 4", "PWRC Fresh", "WP 27 HM-MARC", "WP 28 HM-MARC", 
"WP 30 IT MARE", "WP29 LM-MARC"), class = "factor"), Season = structure(c(2L, 
3L, 2L, 2L, 2L, 3L, 3L, 2L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 
3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 2L, 3L, 3L), .Label = c("", 
"Fall", "Spring", "Spring?"), class = "factor"), Refuge = structure(c(3L, 
2L, 5L, 5L, 5L, 5L, 4L, 7L, 6L, 5L, 5L, 2L, 7L, 7L, 5L, 5L, 5L, 
4L, 6L, 6L, 3L, 5L, 5L, 6L, 6L, 3L, 2L, 6L, 2L, 5L), .Label = c("", 
"Barren Island", "Bishop's Head", "Bishops Head", "Blackwater", 
"Eastern Neck", "Martin", "PWRC"), class = "factor"), Habitat.Type = structure(c(3L, 
3L, 2L, 3L, 4L, 3L, 4L, 3L, 2L, 3L, 2L, 4L, 3L, 3L, 2L, 2L, 5L, 
3L, 4L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 4L, 2L, 3L, 2L), .Label = c("Fresh", 
"High", "Low", "Mud", "SAV"), class = "factor"), Longitude = c(-76.03896, 
-76.26205, -76.05714, -76.06332, -76.14641, -76.23522, -76.03869, 
-75.99733, -76.21661, -76.23491, -76.22003, -76.26163, -75.99354, 
-76.01407, -76.05714, -76.01762, -76.10363, -76.04883, -76.21547, 
-76.23986, -76.03896, -76.01762, -76.05714, -76.2181, -76.23986, 
-76.04883, -76.26163, -76.21661, -76.26205, -76.0235), Latitude = c(38.22447, 
38.33905, 38.40959, 38.39708, 38.41795, 38.43055, 38.23255, 37.99369, 
39.03264, 38.43141, 38.41026, 38.33606, 37.98833, 38.01108, 38.40959, 
38.41913, 38.40351, 38.22694, 39.04036, 39.02677, 38.22447, 38.41913, 
38.40959, 39.03887, 39.02677, 38.22694, 38.33606, 39.03264, 38.33905, 
38.39138), Prey = structure(c(11L, 41L, 35L, 30L, 41L, 41L, 41L, 
3L, 18L, 31L, 40L, 9L, 41L, 38L, 30L, 13L, 35L, 41L, 20L, 27L, 
4L, 40L, 13L, 35L, 41L, 5L, 5L, 15L, 22L, 20L), .Label = c("Hydrobia", 
"Hydrobia genus", "Hydrobia sp.", "Hydrobia spp", "Melampus bidentatus", 
"Ruppia (maritima or rostellata)", "Ruppia genus", "Ruppia maritima", 
"Schoenoplectus pungens", "Schoenoplectus robustus", "Schoenoplectus spp", 
"Schoenoplectus spp.", "Scirpus acutus", "Scirpus acutus?", "Scirpus americanus", 
"Scirpus fluviatilis", "Scirpus genus", "Scirpus genus 1", "Scirpus genus 1?", 
"Scirpus genus 2", "Scirpus genus 3", "Scirpus genus?", "Scirpus heterochaetus", 
"Scirpus meterochaetus", "Scirpus mevadensis", "Scirpus olney?", 
"Scirpus olneyi", "Scirpus paludosis", "Scirpus paludosus", "Scirpus robustus", 
"Scirpus robustus?", "Scirpus species", "Scirpus subterminalis", 
"Scirpus subtermiralis", "Scirpus validus", "Spartina alterniflora", 
"Spartina genus", "Spartina genus?", "Spartina patens", "Spartina pectinata", 
"Zannichallia palustris"), class = "factor"), Density = c(2.36e-05, 
0.000101477, 0.000335244, 1.17e-05, 1.91e-06, 2.8e-06, 1.72e-05, 
1.34e-05, 2.71e-05, 0.000107843, 2.16e-06, 4.46e-06, 1.22e-05, 
6.61e-05, 0.000263052, 3.91e-05, 0.00034925, 3.69e-06, 8.02e-06, 
2.04e-05, 2.9e-05, 2.05e-05, 0.000564046, 0.001912535, 2.04e-05, 
0.001117905, 0.00255132, 9.03e-05, 4.23e-05, 0.000248282), Intakerate_kcals = c(-3.5399430250046e-07, 
7.6382794280604e-14, -5.02872205332896e-06, -1.7549698484651e-07, 
2.70599529637464e-17, 5.81535679492809e-17, 2.19440708445348e-15, 
4.34155540862746e-08, -4.06493587341127e-07, -1.61763139817e-06, 
-3.23994151550826e-08, -6.68988064422799e-08, 1.10402768540446e-15, 
-9.91487886840506e-07, -3.94580269988612e-06, -5.8649138992111e-07, 
-5.23882134070119e-06, 1.00998060784975e-16, -1.2029789281118e-07, 
-3.05994985702607e-07, 9.3958523768985e-08, -3.07494963928282e-07, 
-8.46097103856411e-06, -2.86925082960488e-05, 3.08688633134856e-15, 
3.62058033172122e-06, 8.25888178764606e-06, -1.35448644277712e-06, 
-6.34490870510011e-07, -3.72424640639279e-06), Intakerate_kjs = c(-1.48111216166192e-06, 
3.19585611270047e-13, -2.10401730711284e-05, -7.34279384597799e-07, 
1.13218843200315e-16, 2.43314528299791e-16, 9.18139924135334e-15, 
1.81650678296973e-07, -1.70076916943527e-06, -6.76816976994329e-06, 
-1.35559153008866e-07, -2.79904606154499e-07, 4.61925183573226e-15, 
-4.14838531854068e-06, -1.65092384963235e-05, -2.45387997542992e-06, 
-2.19192284894938e-05, 4.22575886324335e-16, -5.03326383521979e-07, 
-1.28028302017971e-06, 3.93122463449433e-07, -1.28655892907593e-06, 
-3.54007028253523e-05, -0.000120049454710668, 1.29155324103624e-14, 
1.51485081079216e-05, 3.45551613995111e-05, -5.66717127657947e-06, 
-2.65470980221389e-06, -1.55822469643474e-05), Intakerate_kjday = c(-0.12796809076759, 
2.76121968137321e-08, -1.81787095334549, -0.0634417388292498, 
9.78210805250721e-12, 2.1022375245102e-11, 7.93272894452929e-10, 
0.0156946186048585, -0.146946456239208, -0.584769868123101, -0.011712310819966, 
-0.0241837579717487, 3.99103358607267e-10, -0.358420491521915, 
-1.42639820608235, -0.212015229877145, -1.89382134149226, 3.65105565784226e-11, 
-0.043487399536299, -0.110616452943527, 0.033965780842031, -0.111158691472161, 
-3.05862072411043, -10.3722728870017, 1.11590200025531e-09, 1.30883110052443, 
2.98556594491776, -0.489643598296466, -0.22936692691128, -1.34630613771962
)), row.names = c(NA, -30L), class = "data.frame")

lat=unique(profit$Latitude)

## for each location I am calculating the weight for    Fall only
nfall=0
latfall<-c(double())
for(i in lat){
  name = paste0("profit.fall.",round(i,5))
  x = subset(profit,Latitude==i & Season=="Fall")
  if(nrow(x)>=1){
    for(j in 1:nrow(x)){
      x$alpha[j]<- 1 # used to be this  x$Density[j]/sum(x$Density)
    }
    nfall= nfall+1
    assign(name, data.frame(x))
    latfall<-c(latfall,round(i,5))
    print(name)
  }
}

View(profit.fall.38.4.959)

profit.fall.all <- data.frame(matrix(ncol=7,nrow=nfall))

names(profit.fall.all)[1]<-'Id'
names(profit.fall.all)[2]<-'Refuge'
names(profit.fall.all)[3]<-'Season'
names(profit.fall.all)[4]<-'HType'
names(profit.fall.all)[5]<-'Lat'
names(profit.fall.all)[6]<-'Long'
names(profit.fall.all)[7]<-'IntakeEnergy'

View(profit.fall.all)

k=0

lat=latfall

for(i in lat){
  df=as.name(paste0('profit.fall.',i))
  d=get(as.character(df))
  x=0
  for(j in 1:nrow(d)){
    x= x + d$Intakerate_kjday[j]*d$alpha[j]
  }
  k=k+1
  new_row <- c(k,as.character(d$Refuge[1]),as.character(d$Season[1]),as.character(d$Habitat.Type[1]),as.numeric(d$Latitude[1]),as.numeric(d$Longtitude[1]),as.numeric(x))
  #names(new_row)<-c("id","Refuge","Season","HType","Lat","Long","Intakerate_kjday")
  #profit.spring.all <- rbind(profit.spring.all, new_row)
  profit.fall.all[k,] <- new_row
}

View(profit.fall.all)

标签： rfor-loop

解决方案

有问题的代码显然计算（非常低效和不准确）

sum(d$Intakerate_kjday * d$alpha)

但是，您的错误表明，其中一个数据框中缺少一列。

看看new_row这里：

for(i in lat){
  df=as.name(paste0('profit.fall.',i))
  d=get(as.character(df))
  x=0
  for(j in 1:nrow(d)){
    x= x + d$Intakerate_kjday[j]*d$alpha[j]
  }
  k=k+1
  new_row <- c(k,as.character(d$Refuge[1]),as.character(d$Season[1]),as.character(d$Habitat.Type[1]),as.numeric(d$Latitude[1]),as.numeric(d$Longtitude[1]),as.numeric(x))
  #names(new_row)<-c("id","Refuge","Season","HType","Lat","Long","Intakerate_kjday")
  #profit.spring.all <- rbind(profit.spring.all, new_row)

  if (length(new_row) != ncol(profit.fall.all)) {
    # Catch the bad df
    browser()
  }

  profit.fall.all[k,] <- new_row
}

推荐阅读