Hi, thanks for the prompt reply.
The model is a student retention model (binary logistic regression, flat, and hierarchical, by Schoolcode, only on intercepts)
formula = "didNotReturnNextFallIR['1'] ~ pell_amount + studentOfColor + gender + isFirstGeneration + isHEOP + USCitizen + HSGPA + NumAPCourses + applicantTypeCode + WaitListYN + UnmetNeed + hasLoans + isDivisionI + isCampusWorkStudy + GPA_EffectiveEndOfAcademicYear + isSelfDevelopment"
or for the hierarchical case:
formula_h= "didNotReturnNextFallIR['1'] ~ pell_amount + studentOfColor + gender + isFirstGeneration + isHEOP + USCitizen + HSGPA + NumAPCourses + applicantTypeCode + WaitListYN + UnmetNeed + hasLoans + isDivisionI + isCampusWorkStudy + GPA_EffectiveEndOfAcademicYear + isSelfDevelopment + (1|SchoolCode)
This is the dataset (unfortunately I cannot share the actual data without authorization)
All numeric data is scaled ((x-mean(x)/stdev(x)), all non-numeric data (even binary indicators) are category.
<class 'pandas.core.frame.DataFrame'>
Index: 2578 entries, 10798 to 13386
Data columns (total 37 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 AGE 2578 non-null float64
1 EFC 2578 non-null float64
2 pellEligable 2578 non-null category
3 pell_amount 2578 non-null float64
4 studentOfColor 2578 non-null category
5 gender 2578 non-null category
6 isFirstGeneration 2578 non-null category
7 distanceFromHome 2578 non-null float64
8 isHEOP 2578 non-null category
9 USCitizen 2578 non-null category
10 HSGPA 2578 non-null float64
11 NumAPCourses 2578 non-null float64
12 NumHonorsCourses 2578 non-null float64
13 NumIBHCourses 2578 non-null float64
14 NumIBSCourses 2578 non-null float64
15 applicantTypeCode 2578 non-null category
16 WaitListYN 2578 non-null category
17 isMeritScholarship 2578 non-null category
18 UnmetNeed 2578 non-null float64
19 hasLoans 2578 non-null category
20 initialHousingAssignment 2578 non-null category
21 initialHousingRelocationCount 2578 non-null float64
22 SchoolCode 2578 non-null category
23 MAJOR 2578 non-null category
24 isDivisionI 2578 non-null category
25 isLSP 2578 non-null category
26 isHonorsProgram 2578 non-null category
27 isCampusWorkStudy 2578 non-null category
28 GPA_EffectiveEndOfAcademicYear 2578 non-null float64
29 isAcademicProbation 2578 non-null category
30 isSelfDevelopment 2578 non-null category
31 isCareerPlanning 2578 non-null category
32 tutoringClassCount 2578 non-null float64
33 percentCreditsFullTimeFall 2578 non-null float64
34 ISDEANSLIST 2578 non-null category
35 didNotReturnNextFallIR 2578 non-null category
36 DepositDaysDifference 2578 non-null float64
dtypes: category(22), float64(15)
memory usage: 381.9 KB
For the flat model:
model = bmb.Model(formula, dfr, family="bernoulli")
print(model)
Formula: didNotReturnNextFallIR['1'] ~ pell_amount + studentOfColor + gender + isFirstGeneration + isHEOP + USCitizen + HSGPA + NumAPCourses + applicantTypeCode + WaitListYN + UnmetNeed + hasLoans + isDivisionI + isCampusWorkStudy + GPA_EffectiveEndOfAcademicYear + isSelfDevelopment
Family: bernoulli
Link: p = logit
Observations: 2578
Priors:
target = p
Common-level effects
Intercept ~ Normal(mu: 0.0, sigma: 17.0322)
pell_amount ~ Normal(mu: 0.0, sigma: 2.5)
studentOfColor ~ Normal(mu: 0.0, sigma: 7.6971)
gender ~ Normal(mu: 0.0, sigma: 5.096)
isFirstGeneration ~ Normal(mu: 0.0, sigma: 6.744)
isHEOP ~ Normal(mu: 0.0, sigma: 26.5866)
USCitizen ~ Normal(mu: 0.0, sigma: 16.5813)
HSGPA ~ Normal(mu: 0.0, sigma: 2.5)
NumAPCourses ~ Normal(mu: 0.0, sigma: 2.5)
applicantTypeCode ~ Normal(mu: [0. 0. 0.], sigma: [ 8.7191 18.8848 5.3966])
WaitListYN ~ Normal(mu: 0.0, sigma: 8.2321)
UnmetNeed ~ Normal(mu: 0.0, sigma: 2.5)
hasLoans ~ Normal(mu: 0.0, sigma: 5.0546)
isDivisionI ~ Normal(mu: 0.0, sigma: 7.3701)
isCampusWorkStudy ~ Normal(mu: 0.0, sigma: 7.7406)
GPA_EffectiveEndOfAcademicYear ~ Normal(mu: 0.0, sigma: 2.5)
isSelfDevelopment ~ Normal(mu: 0.0, sigma: 8.2321)
start_time = time.time()
idata = model.fit(draws=2000, cores=4)
end_time = time.time()
print(f"Time taken: {end_time - start_time} seconds")
Modeling the probability that didNotReturnNextFallIR==1
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Intercept, pell_amount, studentOfColor, gender, isFirstGeneration, isHEOP, USCitizen, HSGPA, NumAPCourses, applicantTypeCode, WaitListYN, UnmetNeed, hasLoans, isDivisionI, isCampusWorkStudy, GPA_EffectiveEndOfAcademicYear, isSelfDevelopment]
Sampling 4 chains, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 100% 0:00:00 / 0:00:17
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 18 seconds.
Time taken: 32.90929651260376 seconds
For the hierarchical model:
model_h= bmb.Model(formula_h, dfr, family="bernoulli")
print(model_h)
Formula: didNotReturnNextFallIR['1'] ~ pell_amount + studentOfColor + gender + isFirstGeneration + isHEOP + USCitizen + HSGPA + NumAPCourses + applicantTypeCode + WaitListYN + UnmetNeed + hasLoans + isDivisionI + isCampusWorkStudy + GPA_EffectiveEndOfAcademicYear + isSelfDevelopment + (1|SchoolCode)
Family: bernoulli
Link: p = logit
Observations: 2578
Priors:
target = p
Common-level effects
Intercept ~ Normal(mu: 0.0, sigma: 17.0322)
pell_amount ~ Normal(mu: 0.0, sigma: 2.5)
studentOfColor ~ Normal(mu: 0.0, sigma: 7.6971)
gender ~ Normal(mu: 0.0, sigma: 5.096)
isFirstGeneration ~ Normal(mu: 0.0, sigma: 6.744)
isHEOP ~ Normal(mu: 0.0, sigma: 26.5866)
USCitizen ~ Normal(mu: 0.0, sigma: 16.5813)
HSGPA ~ Normal(mu: 0.0, sigma: 2.5)
NumAPCourses ~ Normal(mu: 0.0, sigma: 2.5)
applicantTypeCode ~ Normal(mu: [0. 0. 0.], sigma: [ 8.7191 18.8848 5.3966])
WaitListYN ~ Normal(mu: 0.0, sigma: 8.2321)
UnmetNeed ~ Normal(mu: 0.0, sigma: 2.5)
hasLoans ~ Normal(mu: 0.0, sigma: 5.0546)
isDivisionI ~ Normal(mu: 0.0, sigma: 7.3701)
isCampusWorkStudy ~ Normal(mu: 0.0, sigma: 7.7406)
GPA_EffectiveEndOfAcademicYear ~ Normal(mu: 0.0, sigma: 2.5)
isSelfDevelopment ~ Normal(mu: 0.0, sigma: 8.2321)
Group-level effects
1|SchoolCode ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 17.0322))
start_time = time.time()
idata_h = model_h.fit(draws=2000, cores=1, target_accept=0.95)
end_time = time.time()
print(f"Time taken: {end_time - start_time} seconds")
Modeling the probability that didNotReturnNextFallIR==1 Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Sequential sampling (4 chains in 1 job) NUTS: [Intercept, pell_amount, studentOfColor, gender, isFirstGeneration, isHEOP, USCitizen, HSGPA, NumAPCourses, applicantTypeCode, WaitListYN, UnmetNeed, hasLoans, isDivisionI, isCampusWorkStudy, GPA_EffectiveEndOfAcademicYear, isSelfDevelopment, 1|SchoolCode_sigma, 1|SchoolCode_offset]
Sampling chain 0, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 100% 0:00:00 / 0:01:34
Sampling chain 1, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 100% 0:00:00 / 0:01:28
Sampling chain 2, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 100% 0:00:00 / 0:01:31
Sampling chain 3, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 100% 0:00:00 / 0:01:27
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 363 seconds.
Time taken: 369.7613961696625 seconds
If for example, for the flat model I do
start_time = time.time()
idata = model.fit(draws=2000, cores=4)
end_time = time.time()
print(f"Time taken: {end_time - start_time} seconds")
Modeling the probability that didNotReturnNextFallIR==1
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Intercept, pell_amount, studentOfColor, gender, isFirstGeneration, isHEOP, USCitizen, HSGPA, NumAPCourses, applicantTypeCode, WaitListYN, UnmetNeed, hasLoans, isDivisionI, isCampusWorkStudy, GPA_EffectiveEndOfAcademicYear, isSelfDevelopment]
Sampling 4 chains, 0 divergences ββββββββββββββββββββββββββββββββββββββββ 0% -:--:-- / 0:01:24
It getβs stuck in 0% the clock runs.