Trong lập trình các lệnh ước lượng tùy chỉnh trong Stata, việc đảm bảo chức năng `predict` hoạt động chính xác sau khi chạy mô hình là một yếu tố quan trọng để người dùng có thể trích xuất các dự đoán hiệu quả. Bài viết này sẽ hướng dẫn cách xây dựng một ado-command riêng biệt để tính toán các giá trị dự đoán và tích hợp nó vào lệnh ước lượng chính, từ đó mở rộng khả năng phân tích dữ liệu.
Một Lệnh Ado Để Tính Toán Dự Đoán
Theo quy ước, một ado-command dùng để tính toán các dự đoán cho lệnh ước lượng `mytest` sẽ được đặt tên là `mytest_p`. Trong trường hợp này, chúng ta sẽ xem xét `mypoisson5_p`, ado-command chịu trách nhiệm tính toán dự đoán sau khi chạy lệnh `mypoisson5`.
Cú pháp của `mypoisson5_p` được định nghĩa như sau:
`mypoisson5_p [type] newvarname [if] [in] [, n xb]`
Lệnh `mypoisson5_p` sẽ tính toán số lượng đếm dự kiến khi tùy chọn `n` được chỉ định, và tính toán các dự đoán tuyến tính khi tùy chọn `xb` được chỉ định. `n` là tùy chọn mặc định nếu người dùng không chỉ định cả `xb` và `n`. Điều quan trọng cần lưu ý là người dùng không được phép chỉ định đồng thời cả `xb` và `n`.
Dưới đây là mã nguồn cho ado-command này.
1program define mypoisson5_p
2 version 14
3 syntax newvarname [if] [in] , [ xb n ]
4 marksample touse, novarlist
5 local nopts : word count `xb' `n'
6 if `nopts' >1 {
7 display "{err}chỉ có thể chỉ định một thống kê"
8 exit 498
9 }
10 if `nopts' == 0 {
11 local n n
12 display "số lượng đếm dự kiến"
13 }
14 if "`xb'" != "" {
15 _predict `typlist' `varlist' if `touse' , xb
16 }
17 else {
18 tempvar xbv
19 quietly _predict double `xbv' if `touse' , xb
20 generate `typlist' `varlist' = exp(`xbv') if `touse'
21 }end
Dòng thứ 5 sử dụng lệnh `syntax` để phân tích cú pháp được người dùng cung cấp. Lệnh này chỉ định rằng `mypoisson5_p` yêu cầu tên của một biến mới, cho phép điều kiện `if` hoặc `in`, và chấp nhận các tùy chọn `xb` và `n`. `syntax newvarname` đảm bảo người dùng phải cung cấp một tên cho biến mới không tồn tại trong tập dữ liệu đang được sử dụng. `syntax` sẽ lưu tên biến mới này vào macro cục bộ `varlist`. Nếu người dùng chỉ định một loại biến (ví dụ `double`) ngoài tên biến, loại biến đó sẽ được lưu vào macro cục bộ `typlist`. Chẳng hạn, nếu người dùng nhập:
1mypoisson5_p double yhatthì macro cục bộ `varlist` sẽ chứa "yhat" và `typlist` sẽ chứa "double". Nếu không có loại biến nào được chỉ định, `typlist` sẽ trống.
Dòng thứ 7 sử dụng `marksample` để tạo một biến nhận diện mẫu, với tên được lưu trong macro cục bộ `touse`. Ở đây, tùy chọn `novarlist` được chỉ định để `marksample` chỉ sử dụng các điều kiện `if` hoặc `in` do người dùng cung cấp để tạo biến nhận diện mẫu, mà không tính đến các quan sát không tồn tại trong biến mới.
Các tùy chọn `xb` và `n` chỉ định thống kê cần tính toán. Lệnh `syntax` ở dòng 5 cho phép người dùng chỉ định:
1. Tùy chọn `xb`
2. Tùy chọn `n`
3. Cả `xb` và `n`
4. Không `xb` cũng không `n`
Trong trường hợp (1), macro cục bộ `xb` sẽ chứa "xb" và `n` sẽ trống. Trong trường hợp (2), `xb` sẽ trống và `n` sẽ chứa "n". Trong trường hợp (3), cả `xb` và `n` đều chứa giá trị tương ứng. Trong trường hợp (4), cả hai macro đều trống. Biểu đồ cú pháp và mô tả ngụ ý rằng các trường hợp (1), (2) và (4) là hợp lệ, trong khi trường hợp (3) sẽ là lỗi. Dòng 9 lưu số lượng tùy chọn được người dùng chỉ định vào macro cục bộ `nopts`. Phần còn lại của mã sử dụng `nopts`, `xb` và `n` để xử lý các trường hợp (1) đến (4).
Các dòng 10–13 xử lý trường hợp (3) bằng cách thoát với thông báo lỗi lịch sự khi `nopts` chứa 2. Các dòng 15–18 xử lý trường hợp (4) bằng cách gán "n" vào macro cục bộ `n` khi `nopts` chứa 0.
Tại thời điểm này, chúng ta đã xử lý các trường hợp (3) và (4), và sẽ sử dụng `xb` và `n` để xử lý các trường hợp (1) và (2), vì hoặc `xb` không trống và `n` trống, hoặc `xb` trống và `n` không trống.
Các dòng 20–22 xử lý trường hợp (1) bằng cách sử dụng `_predict` để tính toán các dự đoán `xb` khi macro cục bộ `xb` không trống. Lưu ý rằng các dự đoán được tính toán với độ chính xác do người dùng chỉ định. Các dòng 23–27 xử lý trường hợp (2) bằng cách sử dụng `_predict` để tính `xb` vào một biến tạm thời, sau đó biến này được dùng để tính `n`. Lưu ý rằng biến tạm thời cho `xb` luôn được tính với độ chính xác kép (double precision) và `n` được tính với độ chính xác do người dùng chỉ định.
Lưu Tên Lệnh Dự Đoán Vào E(Predict)
Để người dùng có thể tính toán thống kê `xb` bằng cách nhập:
1predict double yhat, xbthay vì phải gõ:
1mypoisson5_p double yhat, xbCú pháp `predict` tiêu chuẩn này hoạt động vì lệnh `predict` sẽ sử dụng ado-command có tên được lưu trong macro `e(predict)`. Trong dòng 50 của `mypoisson5` (khối mã bên dưới), chúng ta lưu "mypoisson5_p" vào `e(predict)`. Việc bổ sung này là điểm khác biệt duy nhất giữa `mypoisson5.ado` trong khối mã này và `mypoisson4.ado` đã được thảo luận trong các bài viết trước.
1program define mypoisson5, eclass sortpreserve
2 version 14
3 syntax varlist(numeric ts fv min=2) [if] [in] [, noCONStant vce(string) ]
4 marksample touse
5 _vce_parse `touse' , optlist(Robust) argoptlist(CLuster) : , vce(`vce')
6 local vce "`r(vce)'"
7 local clustervar "`r(cluster)'"
8 if "`vce'" == "robust" | "`vce'" == "cluster" {
9 local vcetype "Robust"
10 }
11 if "`clustervar'" != "" {
12 capture confirm numeric variable `clustervar'
13 if _rc {
14 display in red "invalid vce() option"
15 display in red "cluster variable {bf:`clustervar'} is " ///
16 "string variable instead of a numeric variable"
17 exit(198)
18 }
19 sort `clustervar'
20 }
21 gettoken depvar indepvars : varlist
22 _fv_check_depvar `depvar'
23 tempname b mo V N rank
24 getcinfo `indepvars' , `constant'
25 local cnames "`r(cnames)'"
26 matrix `mo' = r(mo)
27 mata: mywork("`depvar'", "`cnames'", "`touse'", "`constant'", ///
28 "`b'", "`V'", "`N'", "`rank'", "`mo'", "`vce'", "`clustervar'")
29 if "`constant'" == "" {
30 local cnames "`cnames' _cons"
31 }
32 matrix colnames `b' = `cnames'
33 matrix colnames `V' = `cnames'
34 matrix rownames `V' = `cnames'
35 ereturn post `b' `V', esample(`touse') buildfvinfo
36 ereturn scalar N = `N'
37 ereturn scalar rank = `rank'
38 ereturn local vce "`vce'"
39 ereturn local vcetype "`vcetype'"
40 ereturn local clustvar "`clustervar'"
41 ereturn local predict "mypoisson5_p"
42 ereturn local cmd "mypoisson5"
43 ereturn displayprogram getcinfo, rclass
syntax varlist(ts fv), [ noCONStant ]
_rmcoll `varlist' , `constant' expand
local cnames `r(varlist)'
local p : word count `cnames'
if "`constant'" == "" {
local p = `p' + 1
local cons _cons
}
tempname b mo
matrix `b' = J(1, `p', 0)
matrix colnames `b' = `cnames' `cons'
_ms_omit_info `b'
matrix `mo' = r(omit)
return local cnames "`cnames'"
return matrix mo = `mo'
end
mata:
void mywork( string scalar depvar, string scalar indepvars,
string scalar touse, string scalar constant,
string scalar bname, string scalar Vname,
string scalar nname, string scalar rname,
string scalar mo,
string scalar vcetype, string scalar clustervar)
{
real vector y, b
real matrix X, V, Ct
real scalar n, p, rank
y = st_data(., depvar, touse)
n = rows(y)
X = st_data(., indepvars, touse)
if (constant == "") {
X = X,J(n, 1, 1)
}
p = cols(X)
Ct = makeCt(mo)
S = optimize_init()
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_evaluator(S, &plleval3())
optimize_init_evaluatortype(S, "gf2")
optimize_init_params(S, J(1, p, .01))
optimize_init_constraints(S, Ct)
b = optimize(S)
if (vcetype == "robust") {
V = optimize_result_V_robust(S)
}
else if (vcetype == "cluster") {
cvar = st_data(., clustervar, touse)
optimize_init_cluster(S, cvar)
V = optimize_result_V_robust(S)
}
else { // vcetype must IID
V = optimize_result_V_oim(S)
}
rank = p - diag0cnt(invsym(V))
st_matrix(bname, b)
st_matrix(Vname, V)
st_numscalar(nname, n)
st_numscalar(rname, rank)
}
real matrix makeCt(string scalar mo)
{
real vector mo_v
real scalar ko, j, p
mo_v = st_matrix(mo)
p = cols(mo_v)
ko = sum(mo_v)
if (ko>0) {
Ct = J(0, p, .)
for(j=1; j<=p; j++) {
if (mo_v[j]==1) {
Ct = Ct \ e(j, p)
}
}
Ct = Ct, J(ko, 1, 0)
}
else {
Ct = J(0,p+1,.)
}
return(Ct)
}
void plleval3(real scalar todo, real vector b, ///
real vector y, real matrix X, ///
val, grad, hess)
{
real vector xb, mu
xb = X*b'
mu = exp(xb)
val = (-mu + y:*xb - lnfactorial(y))
if (todo>=1) {
grad = (y - mu):*X
}
if (todo==2) {
hess = -quadcross(X, mu, X)
}
}
end
end
Ví Dụ Thực Hành
Ví dụ sau minh họa rằng việc triển khai của chúng ta hoạt động bằng cách so sánh các dự đoán thu được sau khi chạy `mypoisson5` với các dự đoán thu được sau khi chạy lệnh `poisson` tích hợp sẵn của Stata.
1clear all
2use accident3
3quietly poisson accidents cvalue kids traffic
4predict double n1
5(option n assumed; predicted number of events)
6quietly mypoisson5 accidents cvalue kids traffic
7predict double n2
8expected counts
9list n1 n2 in 1/5
10 +-----------------------+
11 | n1 n2 |
12 |-----------------------|
13 1. | .15572052 .15572052 |
14 2. | .47362502 .47362483 |
15 3. | .46432954 .46432946 |
16 4. | .84841301 .84841286 |
17 5. | .40848207 .40848209 |
18 +-----------------------+Các kết quả gần như giống hệt nhau, cho thấy rằng `predict` đã hoạt động chính xác sau lệnh `mypoisson5` tùy chỉnh của chúng ta.
✨ Khả năng tích hợp chức năng `predict` vào các lệnh ước lượng tùy chỉnh trong Stata là vô cùng giá trị. Nó không chỉ giúp người dùng dễ dàng trích xuất các dự đoán phù hợp với mô hình độc đáo của họ mà còn biến một lệnh ước lượng cơ bản thành một giải pháp phân tích hoàn chỉnh, nâng cao đáng kể hiệu quả và tính linh hoạt trong công việc khoa học dữ liệu.
1. Hãy suy nghĩ về cách bạn có thể mở rộng `mypoisson5_p` để tính toán khoảng tin cậy cho các dự đoán số lượng đếm. Bạn sẽ cần những thông tin nào từ kết quả ước lượng chính và các hàm nào của Stata sẽ hữu ích?
2. Nếu bạn cần hỗ trợ thêm các loại dự đoán khác (ví dụ: xác suất sự kiện, phần dư) cho lệnh ước lượng tùy chỉnh, bạn sẽ cấu trúc lệnh `_p` của mình như thế nào để xử lý đa dạng các tùy chọn này một cách hiệu quả?
