I am running logistic regressions storing each of the estimates in order to run multiple likelihood-ration test and generate a matrix for p-value of each likelihood test.
Here is a part of my data:
ID pheno cov1 cov2 a_15 a_11 a_3 a_9_f a_9_y ... dpb1_229
11 0 0.0071 -0.0285 1 0 0 0 2 ... 0
22 1 0.0066 -0.0018 1 0 0 0 2 ... 0
23 0 0 0.0085 0 0 0 0.997 0.003 ... 0
27 1 0.0165 0.0084 1 1 0 0.999 0.001 ... 1
35 0 0.0164 0.0123 0 0 0 0 0 ... 0
There are 1,129 columns from a_15 to dpb1_229 indicating different positions.

Using this data, I successfully could run following command for the purpose I mentioned above thanks to the advice I got in a different topic of this forum.

Code:
logit pheno cov1 cov2
est sto E1

local j=2
foreach var of varlist a_15_-dpb1_229{
logit pheno cov1 cov2 `var'
est sto E`j'
local ++j
}

matrix p = J(1130, 1, .)
forvalues i = 2/1130 {
 lrtest E1 E`i', force
 matrix p[`i', 1] = r(p)
}
matrix list p
However, I found that I had to include a_9_f and a_9_y in the same model (logit pheno cov1 cov2 a_9_f a_9_y) in stead of running separately because these two are data of the same position.
They have the same position names, but have additional alphabet characters coding different amino acids.
There are also several similar data in the 1,129 columns, which I also have to deal with in the same manner.
So my question is if there is a way to modify the command above to accomplish this task.

I can also provide additional information about my data if necessary.

I guess this might be a complicated task, but I will really appreciate any suggestions or comments.

Thanks in advance.