Stata：提取值并将它们保存为标量（以及更多）

Question

这个问题是Stata: replace, if, forvalues 的后续问题。考虑这些数据：

set seed 123456
set obs 5000
g firmid = "firm" + string(_n)    /* Observation (firm) id */
g nw = floor(100*runiform())      /* Number of workers in a firm */
g double lat = 39+runiform()      /* Latitude in decimal degree of a firm */
g double lon = -76+runiform()     /* Longitude in decimal degree of a firm */

前 10 个观察结果是：

     +--------------------------------------+
     | firmid   nw         lat          lon |
     |--------------------------------------|
  1. |  firm1   81   39.915526   -75.505018 |
  2. |  firm2   35   39.548523   -75.201567 |
  3. |  firm3   10   39.657866    -75.17988 |
  4. |  firm4   83   39.957938   -75.898837 |
  5. |  firm5   56   39.575881   -75.169157 |
  6. |  firm6   73   39.886184   -75.857255 |
  7. |  firm7   27    39.33288   -75.724665 |
  8. |  firm8   75   39.165549    -75.96502 |
  9. |  firm9   64   39.688819   -75.232764 |
 10. | firm10   76   39.012228   -75.166272 |
     +--------------------------------------+

我需要计算公司 1 与所有其他公司之间的距离。所以，vincenty 命令看起来像：

. scalar theLat = 39.915526
. scalar theLon = -75.505018
. vincenty lat lon theLat theLon, hav(distance_km) inkm

vincenty 命令创建 distance_km 变量，该变量具有每个观测值和公司 1 之间的距离。在这里，我手动复制并粘贴两个数字，即 39.915526 和 -75.505018。

问题 1：提取这些数字的语法是什么？

现在，我可以在 distances_km

. egen near_nw_sum = sum(nw)

将在公司 1 的 2 公里范围内创建工人总和。（或者，collapse 命令可以完成这项工作。）

问题 2：我必须对所有公司都这样做，最终数据应如下所示：

     +-----------------------------------------------------------------+
     | firmid   nw         lat          lon            near_nw_sum     |
     |-----------------------------------------------------------------|
  1. |  firm1   81   39.915526   -75.505018  (# workers near firm1)    |
  2. |  firm2   35   39.548523   -75.201567  (# workers near firm2)    |
  3. |  firm3   10   39.657866    -75.17988  (# workers near firm3)    |
  4. |  firm4   83   39.957938   -75.898837  (# workers near firm4)    |
  5. |  firm5   56   39.575881   -75.169157  (# workers near firm5)    |
  6. |  firm6   73   39.886184   -75.857255  (# workers near firm6)    |
  7. |  firm7   27    39.33288   -75.724665  (# workers near firm7)    |
  8. |  firm8   75   39.165549    -75.96502  (# workers near firm8)    |
  9. |  firm9   64   39.688819   -75.232764  (# workers near firm9)    |
 10. | firm10   76   39.012228   -75.166272  (# workers near firm10)   |
     +-----------------------------------------------------------------+

创建 near_nw_sum 变量是我的最终目标。我需要你的帮助来解决我薄弱的数据管理技能。

Answer 1

以下策略与here 基本相同，并且基于您的“最终目标”。同样，根据原始数据集的大小，它可能很有用。joinby 创建观察结果，因此您可能会超过 Stata 限制。但是，我相信它可以满足您的需求。

clear all
set more off

set seed 123456
set obs 10
g firmid = _n   /* Observation (firm) id */
g nw = floor(100*runiform())      /* Number of workers in a firm */
g double lat = 39+runiform()      /* Latitude in decimal degree of a firm */
g double lon = -76+runiform()     /* Longitude in decimal degree of a firm */
gen dum = 1
list

* joinby procedure
tempfile main
save "`main'"

rename (firmid lat lon nw) =0
joinby dum using "`main'"
drop dum

* Pretty print
sort firmid0 firmid
order firmid0 firmid
list, sepby(firmid0)

* Uncomment if you do not want to include workers in the "base" firm.
*drop if firmid0 == firmid

* Compute distance
vincenty lat0 lon0 lat lon, hav(distance_km) inkm
keep if distance_km <= 40 // an arbitrary distance
list, sepby(firmid0)

* Compute workers of nearby-firms
collapse (sum) nw_sum=nw (mean) nw0 lat0 lon0, by(firmid0)
list

它所做的是形成公司的成对组合来计算距离和附近公司的工人总数。此处无需提取问题 1 中所要求的标量。此外，无需将变量 firmid 转换为字符串复杂化。

以下解决了Stata对观察次数限制的问题。

clear all
set more off

* Create empty database
gen x = .
tempfile results
save "`results'", replace

* Create input for exercise
set seed 123456
set obs 500
g firmid = _n   /* Observation (firm) id */
g nw = floor(100*runiform())      /* Number of workers in a firm */
g double lat = 39+runiform()      /* Latitude in decimal degree of a firm */
g double lon = -76+runiform()     /* Longitude in decimal degree of a firm */
gen dum = 1
*list

* Save number of firms
local size = _N
display "`size'"

* joinby procedure
tempfile main
save "`main'"

timer clear 1
timer clear 2
timer clear 3
timer clear 4

quietly {
    timer on 1
    forvalues i=1/`size'{
        timer on 2
        use "`main'" in `i', clear // assumed sorted on firmid
        rename (firmid lat lon nw) =0

        joinby dum using "`main'", unmatched(using)
        drop _merge dum
        order firmid0 firmid
        timer off 2

        timer on 3
        vincenty lat0 lon0 lat lon, hav(dist) inkm
        timer off 3
        keep if dist <= 40 // an arbitrary distance

        timer on 4
        collapse (sum) nw_sum=nw (mean) nw0 lat0 lon0, by(firmid0)

        append using "`results'"
        save "`results'", replace
        timer off 4
    }
    timer off 1
}

use "`results'", clear
sort firmid0
drop x
list

timer list

尽管效率低下，一些使用timer 的测试表明大部分计算时间都进入了vincenty 命令，您将无法逃脱。以下是使用 Intel Core i5 处理器和传统硬盘驱动器（非 SSD）进行 10,000 次观察的时间（以秒为单位）。计时器 1 是总数，而 2、3、4 是组件（大约）。定时器3对应vincenty：

. timer list
   1:   1953.99 /        1 =    1953.9940
   2:    169.19 /    10000 =       0.0169
   3:   1669.95 /    10000 =       0.1670
   4:     94.47 /    10000 =       0.0094

当然，请注意，在这两个代码中都会重复计算距离（例如，计算公司 1-公司 2 和公司 2-公司 1 之间的距离），这可能是您可以避免的。就目前而言，对于 110,000 次观测，​​这将需要很长时间。从好的方面来说，我注意到与第一个设置中相同数量的观察相比，第二个设置需要很少的 RAM。事实上，我的 4GB 机器被后者冻结了。

另外请注意，即使我使用与您相同的种子，数据也不同，因为我创建了不同数量的观察（不是 5000），这会在变量创建过程中产生差异。

（顺便说一句，如果你想将值保存为标量，你可以使用subscripting:scalar latitude = lat[1]）。