我运行了一个来自 GpGp R 包的模型,该模型输出了整个湖泊一年或一个季节的一系列预测(按纬度、经度坐标)。该模型运行良好,我什至可以绘制预测图,但无法将预测放入数据框或表格中(因此我可以获取数字并将它们与我的实际值的年平均值进行实际比较)。我比 base R 更熟悉 tidyverse,所以这可能是问题的一部分。
这里是我在模型中使用的真实数据:
> head(MI_LatLongDepth)
# A tibble: 6 × 10
year season lake station AvgDens AvgBmss_mg basin depth_m lat long
<dbl> <chr> <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 2015 Summer Michigan MI49b 5.73 2.06 Gree… 40 45.5 -87.0
2 2021 Summer Michigan CS9564 167. 409. Midd… 134. 43.6 -87.3
3 2007 Spring Michigan MI23 37.5 40.8 Midd… 93.4 43.1 -87.0
4 2009 Spring Michigan MI23 24.5 59.7 Midd… 93.4 43.1 -87.0
5 2011 Spring Michigan MI23 14 32.1 Midd… 93.4 43.1 -87
6 2013 Spring Michigan MI23 40.1 35.5 Midd… 90 43.1 -87.0
这里是模型预测平均密度的位置和深度的小标题:
> head(mi_matrix_clean)
lat long depth_m
1 46.03333 -85.60000 0.043701
2 46.03333 -85.59917 3.043701
3 46.03333 -85.59833 0.243698
4 46.03333 -85.59750 0.843704
5 46.03333 -85.59667 1.343704
6 46.03333 -85.59583 1.643707
我设置了各种季节对象,这样我就可以按年份(例如 2022)或年份和季节(2022 春季或 2022.25)获得输出
MI_LatLongDepth$season_numeric = NA
MI_LatLongDepth$season_numeric[MI_LatLongDepth$season == "Winter"] = 0 MI_LatLongDepth$season_numeric[MI_LatLongDepth$season == "Spring"] = 0.25 MI_LatLongDepth$season_numeric[MI_LatLongDepth$season == "Summer"] = 0.50 MI_LatLongDepth$season_numeric[MI_LatLongDepth$season == "Fall"] = 0.75
我创建了一年+季节对象以在模型中使用(year_fraction):
MI_LatLongDepth$year_fraction = MI_LatLongDepth$year + MI_LatLongDepth$season_numeric
我设置模型变量,Y是AvgDens,我预测的,X加上主要变量,深度,和locs_time(坐标和year_fraction)
locs_time <- as.matrix(MI_LatLongDepth[,c("long", "lat", "year_fraction")])
X <- model.matrix(~ depth_m, data = MI_LatLongDepth) #adds intercept
Y <- MI_LatLongDepth$AvgDens
完整模型:
mod_time = fit_model(Y, locs_time, X, covfun_name = "exponential_spheretime")
根据 mi_matrix_clean 的一个子集在新位置设置预测(只有大量的纬度、经度和深度横跨湖面):
rows = sample(1:nrow(mi_matrix_clean), 50000)
mi_matrix_clean_subset = mi_matrix_clean[rows,]
locs_pred = as.matrix(mi_matrix_clean_subset[,c("long", "lat")])
X_pred = model.matrix(~mi_matrix_clean_subset$depth_m)
pred = predictions(mod, locs_pred = locs_pred, X_pred = X_pred)
在新位置对 AvgDens 进行预测:
locstime_pred2022.25 = cbind(locs_pred, rep(2022.25, nrow(locs_pred)))
pred_time2022.25 = predictions(mod_time, locs_pred = locstime_pred2022.25, X_pred = X_pred)
我很容易绘制这些预测,如下所示,但我不确定如何最好地将它们放入数据框或 tibble 中。我可以用来提取 AvgDens 的实际预测值的任何东西。
quilt.plot(locs_pred[,1], locs_pred[,2], pred_time2022.25, zlim = c(-50, 300))