Hierarchical / multilevel proportions test I have data from 16 individuals each responding to 12 trials. Each response can be correct or incorrect. The probability of giving a correct response at random is specified at p=x (I'll have to look that up but it's the same in all trials.)
Ordinarily, I think I would use some sort of a proportions test to see if response success rate is higher than p. I have three questions.

*

*Is a two-proportions test appropriate for these type of data?

*Is there a straightforward hierarchical / multilevel implementation of that to reflect that trial responses are nested under participants?

*Is there a free tool out there for this (in R / python)?

 A: *

*Since the data are nested, the answer is technically no.

2, and a bit of 3.  Straight forward depends on your experience, but here is one approach.
This sounds like exactly the right setup for a mixed effect model. Let's work bottom up.
For subject $j$, we assume that the logit of the success rate is
$$ \operatorname{logit}(p_j) = \beta_{0,j} $$
Here, we assume no other cpvariates are used to adjust this estimate.  Note that the intercept is also indexed by $j$, meaning each subject gets their own estimate.  The next level up the hierarchy is to assume that each intercept comes from some population level distribution
$$ \beta_{0,j} \sim \mathcal{N}(B_0, \sigma_B) $$
The population level mean $B_0$ requires estimation.  Let's do this in R.
library(lme4)
library(tidyverse)

#Simulate data
set.seed(2)
n_subs<-16
n_trials<-12
B0<- qlogis(0.2) #Probabilty of success on the log odds scale
sigma_B<-0.2
b0<-rnorm(n_subs, B0, sigma_B)
p<-plogis(b0)
y<-rbinom(n_subs, n_trials, p)


#Model data
d<-tibble(y=y, n_trials=n_trials, p, subject=1:n_subs)

model = glmer(cbind(y, n_trials-y) ~ 1 + (1|subject), data = d, family = binomial())

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(y, n_trials - y) ~ 1 + (1 | subject)
   Data: d

     AIC      BIC   logLik deviance df.resid 
    65.6     67.1    -30.8     61.6       14 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-0.99949 -0.58779  0.04817  0.50596  1.71082 

Random effects:
 Groups  Name        Variance Std.Dev.
 subject (Intercept) 0.1829   0.4277  
Number of obs: 16, groups:  subject, 16

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -1.144      0.207  -5.525  3.3e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Note that the intercept estimate for this model captures our population level mean
> confint(model)
Computing profile confidence intervals ...
                2.5 %     97.5 %
.sig01       0.000000  1.0026083
(Intercept) -1.621016 -0.7514589

These intervals are on the log odds scale.  If you wanted to test if the population level mean success rate was $p$, you would look to this confidence interval and determine if $\operatorname{logit}(p)$ was in this interval.
I think this is the easiest since its just a model call in R, but perhaps you could also weight each observation by their observed variance and calculate a weighted mean.  To obtain a confidence interval for that would be an exercise in the manipulation of random variables and their variances.
EDIT:
Libraries like rstanarm and brms make Bayesian analysis easy, but I'm fond of writing my own little Stan program. Here is a stan program for a hierarchical beta binomial model.
model_code<-cmdstanr::write_stan_file('
data{
  int n;
  int n_subject;
  int n_trials[n];
  int y[n];
  int subject[n];
}
parameters{
  real<lower=0, upper=1> mu;
  real<lower=0> kappa;
  vector<lower=0, upper=1>[n] p;
}
model{
  mu ~ uniform(0,1);
  kappa ~ cauchy(0,1);
  p ~ beta_proportion(mu, kappa);
  for( i in 1:n){
    target += binomial_lpmf(y[i] | n_trials[i], p[i] );
  }
}')

Fitting this model to the data I've simulated above results in the following posterior credible interval for the population probability
> fit$draws('mu') %>% 
+   posterior::as_draws_df() %>% 
+   tidybayes::spread_draws(mu) %>% 
+   tidybayes::mean_qi()
# A tibble: 1 x 6
     mu .lower .upper .width .point .interval
  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 0.263  0.187  0.350   0.95 mean   qi   

