Hi. I've been getting the conformability error r(503) on running the code below. Obviously my matrix specifications are off, but I can't seem to figure out where. Any help would be much appreciated. Thank you.

Code:
local time "r30 r90 r365"
local total "total"
local base_char "ib3.stratumn i.dschrgyr"  
local pati_char "ib1.racecat1n_mice c.age##c.age ib1.stkcat1n ib1.sscat2n_mice ib2.sexn ib0.hisstktian"
local hosp_char "ib2.hrucacat2n ib2.medschcatn"
local other_char "ib2.rucacat2n_mice ib0.dc_refneedn_mice ib0.dc_refhhn_mice ib0.pcpfln ib2.dualcatn ib0.medhisdpn ib0.medhishtn ib0.medhissmn ib0.anycvdn"


*******************************
* coefplots for total payments
*******************************

cap program drop multireg
program multireg
syntax varlist(numeric) [if]
local timestamp=subinstr("`c(current_time)'", ":", "_", .)
local indvar i.tx
*four specifications
local spec1 ib3.stratumn i.dschrgyr
local spec2 ib3.stratumn i.dschrgyr ib1.racecat1n_mice c.age##c.age ib1.stkcat1n ib1.sscat2n_mice ib2.sexn ib0.hisstktian

*placeholder for results
mat m = J(9,2,0)
loc b = 6   
foreach i in `*' {
loc n = 1
while `n' <= 2 {
display " "
display "This is a regression of `i' on `indvar' with `spec`n''."


**********************************************
**     30-day model -   TPM (gamma)    **
**********************************************
if "`i'" == "r30_total_w"   {
tpm `i' `indvar' `spec`n'', first(logit) second(glm, fam(gamma) link(log)) cluster(hospclus)
margins, at(tx=(0 1))
mat a = r(table)  
sca a1 = a[1,1]  
sca a2 = a[1,2]  
mat m[`b', 1] = a1, a2
loc b = `b' + 1
margins, dydx(tx) post
est store m`i'_`n'
outreg2 `*' using "./table/`i'_`timestamp'.doc", ctitle("`i'_`n'") alpha(.01, .05, .1)
}


******************************************************************
** 90-day model - TPM(inverse guassian) **
******************************************************************
else if "`i'" == "r90_total_w" {
tpm `i' `indvar' `spec`n'', first(logit) second(glm, fam(igaussian) link(log)) cluster(hospclus)
margins, at(tx=(0 1))
mat a = r(table)
sca a1 = a[1,1]
sca a2 = a[1,2]
mat m[`b', 1] = a1, a2
loc b = `b' + 1
margins, dydx(tx) post
est store m`i'_`n'
outreg2 `*' using "./table/`i'_`timestamp'.doc", ctitle("`i'_`n'") alpha(.01, .05, .1)
}


**********************************************
** 365-day model - GLM (gamma)    **
**********************************************
else if  "`i'" == "r365_total_w" {
glm `i' `indvar' `spec`n'', fam(gamma) link(log) cluster(hospclus)
margins, at(tx=(0 1))
mat a = r(table)
sca a1 = a[1,1]
sca a2 = a[1,2]
mat m[`b', 1] = a1, a2
loc b = `b' + 1
margins, dydx(tx) post
est store m`i'_`n'
outreg2 `*' using "./table/`i'_`timestamp'.doc", ctitle("`i'_`n'") alpha(.01, .05, .1)
loc n = `n' + 1
}
loc b = `b' + 1
}
putexcel set "$output/table/coefplot.xlsx", sheet(Total Payments) replace
*rownames
matrix rownames m = "0" "0" "30-day" "Base Covariates" "+Patient Covariates" "90-day" "Base Covariates" "+Patient Covariates" "365-day" "Base Covariates" "+Patient Covariates" 
*main table
putexcel A1 = mat(m), nformat(number) rownames
mat list m


**********************************************
** format table     **
**********************************************
*title
putexcel A1 = "Post-acute total payments for COMPASS patients versus control patients by timewindows"
*subtitles
putexcel B2 = "Adjusted mean spending", bold
putexcel A4 = "Spending", bold border(bottom)
putexcel B3 = "Control", bold
putexcel C3 = "Treatment", bold
putexcel B4 = "N=201", bold border(bottom)
putexcel C4 = "N=140", bold border(bottom)
putexcel D4 = "Differences in Spending", bold border(bottom)
putexcel (A1:I1), merge
putexcel (B2:C2), merge hcenter
putexcel (D4:I4), merge hcenter
end

multireg r30_total_w r90_total_w r365_total_w

*generate coefplots
coefplot ///
(mr30_total_w_?, label("30-day")) (mr90_total_w_?, label("90-day")) (mr365_total_w_?, label("365-day")), ///
keep(1.tx) aseq swapnames xline(0) ///
mlabel(cond(@pval<.001, "$"+string(@b, "%9.0f")+"***", cond(@pval<.01, "$"+string(@b, "%9.0f")+"**", cond(@pval<.05, "$"+string(@b, "%9.0f")+"*", cond(@pval<.1, "$"+string(@b, "%9.0f")+"+", "$"+string(@b, "%9.0f")))))) mlabposition(1) mlabsize(medium) ///
yscale(off) groups(mr30_total_w_? = "{bf:30-day}" mr90_total_w_? = "{bf:90-day}" mr365_total_w_? = "{bf:365-day}", nolab) ///
coeflabels(*1 = "base" *2 = "patient" *3 = "hospital" *4 = "other") legend(off) grid(none) ///
note("+ p < .1, * p < .05, ** p < .01, *** p < .001", size(vsmall))

graph export "./graph/totpmt_adjusted.png", replace

*put into excel (to combine with other statistics)
putexcel set "$output/table/coefplot.xlsx", sheet(Total Payments) modify
putexcel D5 = picture("$output/graph/totpmt_adjusted.png")


log close