Hydroxychloroquine statistics.ipynb

Kernel: SageMath 9.1

This notebook is the data analysis I did for this post.

Boulware et al.

The incidence of new illness compatible with Covid-19 did not differ significantly between participants receiving hydroxychloroquine (49 of 414 [11.8%]) and those receiving placebo (58 of 407 [14.3%])

n( (49/414) / (58/407))

0.830543061802432

The study used Fisher's exact test to calculate the p-value. In this case, p is the answer to the question, "Assume that there are $821$ patients, $107$ of whom are sick. We pick $414$ of them at random. We would expect $\frac{107}{821} \cdot 414 \approx 53.96 \approx 54$ to be sick. What's the chance that 49 or less are sick, or 59 or more are sick?"

($49$ is the number we actually observed, and $59$ goes the same distance in the opposite direction, since it's a two-tailed test.)

For instance, the chance that the first 49 patients we choose are sick is $\frac{107}{821} \cdot \frac{106}{820} \cdots \frac{59}{773}$. Of the $772$ patients not yet chosen, $58$ of them are sick. The chance that the next $414 - 49 = 365$ patients we choose are healthy is $\frac{714}{772} \cdot \frac{713}{771} \cdots \frac{350}{408}$. But it was arbitrary to make the first $49$ patients sick, so we multiply by ${414 \choose 49}$.

Let's write a function that does just that.

def chance(sick_patients):
    prob = 1
    for i in range(sick_patients):
        prob = prob * (107-i)/(821-i)
    for i in range(414-sick_patients):
        prob = prob * (821-107 - i)/(821- sick_patients -i)
    prob = prob * binomial(414,sick_patients)
    return prob

Let's check that the chances sum to 1

sum = 0
for i in range(822): #0 to 821
    sum += chance(i)
print(n(sum))

1.00000000000000

sum = 0
for i in range(50): #0 to 49
    sum += chance(i)
for i in range(59,822): #59 to 821
    sum += chance(i)
print(n(sum))

0.350871118206300

And we get p=0.35, just as the paper says.

Now let's do the Bayesian analysis I describe in the blog.

ineffective_proportion = 58/407
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5

print(n(very_effective_proportion))
print(n(effective_proportion))
print(n(ineffective_proportion))
print(n(dangerous_proportion))
print(n(very_dangerous_proportion))

0.0712530712530713
0.128255528255528
0.142506142506143
0.156756756756757
0.213759213759214

#Motorcycle versus stock market
print( 10000 * 0.97^50)
print( 10000 * 1.07^50)

2180.65375347407
294570.250630714

f(x) = x^49 * (1-x)^(414-49) * binomial(414,49, hold=True)

show(f)

latex(f)

x \ {\mapsto}\ -{\left(x - 1\right)}^{365} x^{49} {414 \choose 49}

print(n(f(very_effective_proportion)))
print(n(f(effective_proportion)))
print(n(f(ineffective_proportion)))
print(n(f(dangerous_proportion)))
print(n(f(very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

0.000171386510778484
0.0503315749025967
0.0214397799798485
0.00504589206399895
1.61492837223403e-7

print(n(100*0.2*f(very_effective_proportion)))
print(n(100*0.2*f(effective_proportion)))
print(n(100*0.2*f(ineffective_proportion)))
print(n(100*0.2*f(dangerous_proportion)))
print(n(100*0.2*f(very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

0.00342773021556968
1.00663149805193
0.428795599596971
0.100917841279979
3.22985674446806e-6

total = f(very_effective_proportion) + f(effective_proportion) + f(ineffective_proportion) + f(dangerous_proportion) + f(very_dangerous_proportion)

n(total * 0.2 * 100)

1.53977589900120

print(n(f(very_effective_proportion)/total))
print(n(f(effective_proportion)/total))
print(n(f(ineffective_proportion)/total))
print(n(f(dangerous_proportion)/total))
print(n(f(very_dangerous_proportion)/total))

0.00222612278695435
0.653751950985142
0.278479225369819
0.0655406032432649
2.09761481950923e-6

print(n(100*f(very_effective_proportion)/total))
print(n(100*f(effective_proportion)/total))
print(n(100*f(ineffective_proportion)/total))
print(n(100*f(dangerous_proportion)/total))
print(n(100*f(very_dangerous_proportion)/total))

0.222612278695435
65.3751950985142
27.8479225369819
6.55406032432649
0.000209761481950923

print(n(100*f(ineffective_proportion)/total) + n(100*f(dangerous_proportion)/total)+ n(100*f(very_dangerous_proportion)/total))

34.4021926227904

Digression: what if the data was different? That is, suppose the hydroxychloroquine results were the same as the control?

n(58/407 * 414)

58.9975429975430

n(59/414)

0.142512077294686

n(58/407)

0.142506142506143

Specifically, suppose 59 out of 414 of the hydroxychloroquine group got infected. That's 14.3%, almost exactly the same as the placebo group's 14.3%. Now follow the same procedure as before:

fake_f(x) = x^59 * (1-x)^(414-59) * binomial(414,59, hold=True)

show(fake_f)

print(n(fake_f(very_effective_proportion)))
print(n(fake_f(effective_proportion)))
print(n(fake_f(ineffective_proportion)))
print(n(fake_f(dangerous_proportion)))
print(n(fake_f(very_dangerous_proportion)))

1.96804363176439e-7
0.0388781457121887
0.0560069722228632
0.0404266781168745
0.0000579218083589014

total = fake_f(very_effective_proportion) + fake_f(effective_proportion) + fake_f(ineffective_proportion) + fake_f(dangerous_proportion) + fake_f(very_dangerous_proportion)

print(n(fake_f(very_effective_proportion)/total))
print(n(fake_f(effective_proportion)/total))
print(n(fake_f(ineffective_proportion)/total))
print(n(fake_f(dangerous_proportion)/total))
print(n(fake_f(very_dangerous_proportion)/total))

1.45382645519118e-6
0.287199307235300
0.413732788128065
0.298638572810091
0.000427878000088801

Another digression: What if our priors were different?

ineffective_proportion = 58/407
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5

f(x) = x^49 * (1-x)^(414-49) * binomial(414,49, hold=True)

show(f)

prior_very_effective_proportion = 40
prior_effective_proportion = 40
prior_ineffective_proportion = 10
prior_dangerous_proportion = 5 
prior_very_dangerous_proportion = 5

posterior_very_effective_proportion = prior_very_effective_proportion * f(very_effective_proportion)
posterior_effective_proportion = prior_effective_proportion * f(effective_proportion)
posterior_ineffective_proportion = prior_ineffective_proportion * f(ineffective_proportion)
posterior_dangerous_proportion = prior_dangerous_proportion * f(dangerous_proportion)
posterior_very_dangerous_proportion = prior_very_dangerous_proportion * f(very_dangerous_proportion)

total = posterior_very_effective_proportion + posterior_effective_proportion + posterior_ineffective_proportion + posterior_dangerous_proportion + posterior_very_dangerous_proportion

print(n(posterior_very_effective_proportion/total))
print(n(posterior_effective_proportion/total))
print(n(posterior_ineffective_proportion/total))
print(n(posterior_dangerous_proportion/total))
print(n(posterior_very_dangerous_proportion/total))

0.00303372982676281
0.890924258370064
0.0948769242524658
0.0111647302255928
3.57325114785736e-7

print(n(100*posterior_very_effective_proportion/total))
print(n(100*posterior_effective_proportion/total))
print(n(100*posterior_ineffective_proportion/total))
print(n(100*posterior_dangerous_proportion/total))
print(n(100*posterior_very_dangerous_proportion/total))

0.303372982676281
89.0924258370064
9.48769242524658
1.11647302255928
0.0000357325114785736

End of digressions.

We can also integrate the likelihood function to get probabilties. (Well, we should multiply by the prior, but we assume the prior is uniform, i.e. $y = 1$.)

plot(f,0,1)

plot(f,0,0.2)

total_area = numerical_integral(f, 0, 1)[0]
print(total_area)
at_least_effective_area = numerical_integral(f, 0, effective_proportion)[0]
print(at_least_effective_area)
not_effective_area = numerical_integral(f, effective_proportion,1)[0]
print(not_effective_area)

#To check that the pieces add to the total
print( (at_least_effective_area + not_effective_area) / total_area ) 

print(at_least_effective_area / total_area)

0.0024096385542168746
0.0016948379490692445
0.0007148006051476249
0.9999999999999979
0.7033577488637344

So the probability of effectiveness is roughly 70%.

But really we should do this in 2 dimensions, since the control group's infection rate is also determined by the data, so it may vary.

var('x y')
f_2d(x,y) = x^49 * (1-x)^(414-49) * binomial(414,49,hold=True) * y^58 * (1-y)^(407-58) * binomial(407,58,hold=True)
show(f_2d)

plot3d(f_2d*100, (0,1),(0,1))

Now we need to do some integrals. Symbolic integration is theoretically possible, but the result is a really long expression. You don't believe me? Here it is:

var('x y')
f_2d.integrate(x,0,y)

y |--> 27101836261352403057669001965570401374968441258953491312131661547410800/83*(14589930354045020136586730464187821267314045147212760459986525400*y^415 - 5338187682074805374853804340731522286150712533029113262503523997500*y^414 + 973902579828020424659768642773459421832648639076701273857422581510000*y^413 - 118128236945984623013404413556597300211754734079077652326400433654852500*y^412 + 10716616649070517405542085070246440232835734829373453393289942990887665000*y^411 - 775626892116703286818382782418309553066376006033643245568693722683411582300*y^410 + 46651397423156236811081213808045024463161001829896390809266175007364119620000*y^409 - 2398414332997007495520961381179642502433611015296049355362521618524608210057500*y^408 + 107592749120047574577940326529675461448238009624754961007884323614305397295405000*y^407 - 4278357650354075681061891489993244354944590739044135489733055374065396803947742500*y^406 + 152685605271994242626657172443867547181795399273699547442977010258250911717380471600*y^405 - 4939777879833702122782996167443082042318324376817316029166071028325838171054540347500*y^404 + 146085044098903379402798382936939582581513126359575117706206537259521933803047175835000*y^403 - 3976638401926118608967220657907785701963925725753309733977657057788429158225734440702500*y^402 + 100233388034583806135250728887884451010100346615068817348908867884804752334445017153040000*y^401 - 2351324933209284216223779223616437394020438981069591851779378677275692282637578434884588840*y^400 + 51564143272133425794381123324921872675886819760298066924986374501659918478894263922907650000*y^399 - 1061241282114669180504487999426351713270753109435599485275092763188094192743793252369023525000*y^398 + 20569012389869876878896390410376073005157065304963374070604736646040274009485913885547052100000*y^397 - 376603745057959113365216353873338309558404061582462712465767931369200173747045456002003767975000*y^396 + 6531739080109028592006137207887255926892948215415533039393627706896497595712048126527917502894000*y^395 - 107579495182723939010833496253908122001128416492114063928374016891143365948230524346457814564525000*y^394 + 1686432391751754776517734067576878466752826863020491761345836311977562716534741192206975058387400000*y^393 - 25213997335375420598914600108173873163244710327007461062295628773099320832538440324790154269694225000*y^392 + 360218386484417134182958967274320346521444889735712730163640364006477765142096490164393354860644350000*y^391 - 4925977198805520094128010468272041653543668573677160412866495730599762759401114577447049271479443264600*y^390 + 64582220344236895835611448298682294172423418318132437289509327059277095303202273380411185821452855140000*y^389 - 812955515369880077390026469892662544369191351377476387724788513615304489760269854233314574786255198565000*y^388 + 9838892405654074646478186726430725057101627038118102933903623553854327427567444607195810870779683625010000*y^387 - 114630919122118881661518548747040265429885393664849461009181407921219182488468693784694520096142179622135000*y^386 + 1287201011814920387530070103748583053263585802316564129441280682402635765543241186935042465661407456993283200*y^385 - 13946299605076390875478057092655703905802333573066100352376442877426070888855807466813432628417340958140851875*y^384 + 145944566102209333234767761951342196486829642195428016220169146351758987265102026963415764946988779530622805000*y^383 - 1476568619396293523038953323140745073594618890684015053205152527094795294050757370846190960302554494094566232500*y^382 + 14456101481111280667430308495754350255794653866922149537416262474623306960992222695707682750456320145055524770000*y^381 - 137073189013866520018819822158502019932971398030933774459631750145655628884036856870555757845323081940834382275900*y^380 + 1259819547225070478625037063813848819964424898349391330785533675481179702584155808528854590574869134812418025315000*y^379 - 11231814611972062280154046655513558073086236593742596484125481161834961983048832541203446857812884834040841984452500*y^378 + 97205451825404451329234909991427682823783607276181912780107054053395639153031303489359776296179862715602383959260000*y^377 - 817197960691764549605615825752401344164840379255348474196591483677615998305138165770920885298469537351114722593672500*y^376 + 6677923815316915428830531030495089810823298288506906304874759846884497240550041053835759895756346800690015882442800024*y^375 - 53076298747342446060305267336741168037960706854817999602622365745340296151821850361042991867438683079046281835966887500*y^374 + 410543442934250494306367360832181014819859990670204641889912408499767376898580654879978653295592356746234664097735450000*y^373 - 3092139501205024087492925310433879500155158278066740750828629473846178657249568562095833218954718151314679050752129862500*y^372 + 22689833252410108141157349018929908236129960743686993022375432056799232068119475658054844786772357206387987957320088425000*y^371 - 162291587499265596915229420541344246351288477773359020181844599343920561334799788323928499729820227751708821250763162199500*y^370 + 1132044536344588586942568834222643971249708717353935901350913194433851447039018200633506207612272591363398399675406999200000*y^369 - 7704330158282716648797220169124748178756617878610612215767653967969651230481606811677303954917218952473230170223857150950000*y^368 + 51180264103251398337241478943667945667053499694966410305290300746348800272245660509180564147379236528963856008190091918300000*y^367 - 332009632291076415253503759274416478808929246181147385329941357116620557314195823462560368704246291878171408728494875860550000*y^366 + 2104049645603108705893656481240652727294549090209686684806558772530403680916746967582788054951140226491057120433030321234398400*y^365 - 13031302952123325978244082844787979278210018403730808370266865700909457186673291456465797059793543823573006486844842770088610000*y^364 + 78905795836823885675428358382490298769987935127824564181395511874928338557322023750170308367565314942736827157258855120261060000*y^363 - 467278099509230327663738906028362621492091066469733094846553563208458401883892995750578560742958590559721440882966426556287990000*y^362 + 2707307770132887561693512209011356721134848333402971922865966227711917807652324331181560734818544939580238966814699911817656240000*y^361 - 15351118720231788860550116627167375274673260466640659633943681532901248696107783080105303295899639338970353992083493626246468874800*y^360 + 85216118641955214725298100656181808627812170119234859488067790443283891488541016302216068833187412964041559168151506600134596420000*y^359 - 463251458251040572610606860917911605061836314980733225417565816998134145188441805322999643817225958419783112813986806755833133170000*y^358 + 2466915811484243103367910504793422448591631383900104721696246365569509484451769884194689110685917274181346718406328661263652303480000*y^357 - 12872381737900474189964750686290911585564378596651474833040702029139814796535399330501510752336023963673653541735632019893116968630000*y^356 + 65834074037555213913521130552196560773945593746705401991179432687623233083351185984390261751947281995622212310916094353008290614773600*y^355 - 330100229990566397447740132006070325914557708758197990210009867289636267578950155712691284208210241644715895061796800781326880907410000*y^354 + 1623141110786109541612720554051232029316505637655224111963831759035067923411144108433564161626188127791179099103094421679762113484960000*y^353 - 7828713477118122766114492141989796474732328936596551036775770252283181368151720900701009612345574312789354494091264260842608040054015000*y^352 + 37046988847558979813436385891296444286411391349535245789414684869066564821937345914713324689760680722345064999104244835497355995811165000*y^351 - 172046216208063902253598576079180687266094501427241681446041796531945127033077034427928679859248601274570481855840113016049721244547050260*y^350 + 784269024131575485657195905544585055263520128933766456187999853948694900515459260872297967579278593984106129845585302886394926225468450000*y^349 - 3510004157537915315826644790679490268562514538272751049479926778152352957783868988303532561968832585384948479821439009328280877790095137500*y^348 + 15426405795614170618430512433150782285599551096670510884695131657972182197576637970970907815155781252532662194922974700092390242004747675000*y^347 - 66592525546209299227813947425990050748580218500115865848445822800492381884060096031259562223428820478544174966346940169926348577983048112500*y^346 + 282407457405406691085540092736240042959279157677344369222957655815123303645461647831779682181146421454682410904902203910690555526406288012000*y^345 - 1176794095047300441583313254922651096142345081584168161548743671003282600250444612187986290287675116424430984680904921512059093051916146512500*y^344 + 4819252008288944665531663805873714012773413191249450566342474081251538267692296983246039093559050476785764984883705869049384857260228028575000*y^343 - 19399583737123714249029347145177536832881328699326224964625983445153342320147727089611619116542137862152230842579009698533771079459891718487500*y^342 + 76774194750103931194621016220417275829391337450642550918102925052183299415837069750216958526367759985177553335313871147435068363286939311050000*y^341 - 298760004676375027308762276885028497953973796335871018437422406111890234679722677191520745550196046709378580196726695993377618382640246069600100*y^340 + 1143368128438541286005432405715875248341166601409597982375789832426345953245639495850176428786264951309762770043364858909603833741978274431865000*y^339 - 4304038996052114359918474533942596401876294209893941594226485271368690516615251942209238258073432857326088078994793613779789333571634270512727500*y^338 + 15938993077368067421893935366054481630687285382633943945384728838777821260343722326044894202004878949385631817761134807113285128478930473590160000*y^337 - 58077285693563973514849289280795034992406704148695465056271307585810431527998209751836050419066828176071011293777848003661461134017242181855285000*y^336 + 208246076343609089510558108310719391142922188129591022762964461707688030899795371802702876905621116433123876913689406035815364519944809877147965200*y^335 - 734911441520887034117710334551856800080228280907544278047387890202613837146372555735309276726701267518758858147944704711668216439037876950569305000*y^334 + 2552946594041795228457830934382587293210654574590420421829012880431375453612990373794349108371346056975001473979398041094631228980494281697180180000*y^333 - 8730843884304860666810768232881428506491467428321477360021205720676890382241058202787257425744203841750151107067226447768796356071011780115421795000*y^332 + 29399242117111144981639021904592535696006227395668608702151681930393223462341275830015616070632834861171635277228528591062645304301266253112007270000*y^331 - 97484952176635461885462640512230172464215159226329160926767561070509950354931351380591175270788969813634045232452858209744303115806241549668481575400*y^330 + 318357587514000023970338162979079237141919866845013588380339387959112158506950479520980582279667095878270508923398745266854230858150021884122700840000*y^329 - 1024052471106587335304923966849336356058801741574524891050593382399595743675059274800614263668100728825425016981626514816015443814846112235965336165000*y^328 + 3244968047170665123104121038128644922534951973996784739559261432468888597833840871330942834942733418724546662095240204579634147512264661458307729110000*y^327 - 10130487376599975052927834132234594624605135963764209172280221317005564910153424977055403289662847257691139177842585044827279294392824869726065202835000*y^326 + 31162418194767841208965292440408201228836639776945247332823425409192913106186822760189297996263396405402322689421006221995872881787651640289656160351520*y^325 - 94462797658301987791726801379103343672068738300925297801217394042416789067475487864330433740877799950238522261803094257593309551274296749114367770025000*y^324 + 282205859634991831612543361126372272103633541237143095894497084853600513738531562479600167780999043221913280802053826080633451632922565897138817914650000*y^323 - 830983453441325877237776699417697875570236493446538471010514207082937277130500722897504280864246130777569331057380095170112106175765676158545509413475000*y^322 + 2412038668595648095446193227341408190294412177438245689831804742929408712142830758741135246742696795851482478005222319062076579415995923002255249435600000*y^321 - 6902159730914606035228940626304488154009908214721395750136793631321453726250558899133226850309011084470931323820799169655562622367661228220039026435930250*y^320 + 19473177924210487246727418694902944635137672079151273276248007110311311453371482787523210549461159799447768625199746873636383574078041082752461203110775000*y^319 - 54172757273371786329797267899511008950457069921028655615253768337642609969696567123132795839324761186228764620923798190157936551747739848021508406280837500*y^318 + 148613244893640543300413389540588625081882229026902277285597989854093021053802235128496407410371900062291902274306630301874074058010345097784878401367350000*y^317 - 402073789622581573468983789711222370486684246982738199948079793608892060125711781710804482380976449602746589841293916446756161430633473939758800761198412500*y^316 + 1072911570175102118714568299300524049956467652837493401105898258145683523997677127729631160984614490406618011292057188598366219213121501126360817675660119400*y^315 - 2824030186049273171460509200482791356187051060118867429441724698830514352201113017507440283886948810830940231210332150414599865952428495507440539051388297500*y^314 + 7332606850301965445761773184339456000203025267168252660245605795723056998686445677482201646717174360165433233302377182884441223078134202956361301923319060000*y^313 - 18783073914301950881555371133981817229297760672393060453658141134494844693641695348663721600576023031941072265752377001427194331670678661101714775112823677500*y^312 + 47471048542255091295506500679452438399447073596466062754261411355990186267531744514629212791166411907092227655567100717433359307695027098475716441281927365000*y^311 - 118380106950027281583642339989303652145054329059869332713553088208205204593604277542149549353971392347483355177574105448073216846848327023360134539333189743300*y^310 + 291305996830072629164119576330764925288251173054450394355265888428324463889738803501778521033464438639774025380487051936090287533002844174232974411377391870000*y^309 - 707413224672519194824422072257781518865882680869059337014391440611944264157048373074731060105267401650071459680354389291699290009520288327103992448019906807500*y^308 + 1695436255440254110795267087401747816840345679592801639417046731716357190462531424719822323567528314776247993486538822254662576345654487622903596811721977155000*y^307 - 4010560206211567101026436946960199706208743592198610557905065961674427833852109698356097065034311222268129472869653610873377184478497134345489922793309716242500*y^306 + 9364269577162125795111981063400125311045250016049787083876341197428553764311143526987279275306342888452941591258336371427030555091134173355836023604670936857600*y^305 - 21583276030221204768694710976942411867785215742006643861085894874051854226508392241922583144765019504191190500833247726639915286338696074346554846368340508360000*y^304 + 49109330561927173132352982166300839073046417336597861043658900830822370814299623271814420677338983558334127938246701333693649774573572491653217203476478686160000*y^303 - 110316830907424695691058057196461818628045047957752312891808230639102506131787927691321036004926716890298772596203887872578722911241736319452733329464437074440000*y^302 + 244668418120261381385798285973058525525113374417928018310840286313407271861443971819674072510314823751800949748251462454238244694210558397611202928995789447840000*y^301 - 535795468330546887683433226189870859766604073931851585837747665833745878156636333766614372528312905776009088534207137387830539851950779637620985289523504016630720*y^300 + 1158593591419415770276135738791382611074135207729468023810145165999623482347702832213587762069642138542286980742288527969302574354062636167154131493318613806000000*y^299 - 2473999067136090274560562780036632443557633234402186194019650024489128890158573497422526887909336512838060230690188814198432454633003996021361562663025988541000000*y^298 + 5217098301022245495951715793409441422949361597764628444981687365458100253902848135408161477031661406288421026334751872525279693180320930262018630572626623784000000*y^297 - 10865350734291701888982566400429344008237060454531438169867465462259844893426613569310905002088509229613680390384143272103735890643174060955355005682980822569000000*y^296 + 22349473984976965105880072175527206590502536555287039598222298449150080953665319030026634424634940340778241562322366567839142164211111932839692940503161908931760000*y^295 - 45406989156805817536051110824934751681806117657091437200245027840904090091813148167065131923328563364528796562569546401601562936930289395232158522372126467871000000*y^294 + 91123923632088125294191307730585781532088727584879744278990021742155648375379287106601698194256639106972465456624209160893238931450546636165560447627339123236000000*y^293 - 180641837584392539058389768323587734138077328079848895002778999833493150926568377956981131243258671560789970153864212831229472064827582969939629644465389192429000000*y^292 + 353754428862196939918580384348132237144803688932336833956223679310330340790240213389649558202350924708226982899696838917692035638758904860542456191012699864034000000*y^291 - 684390395876820237060136853371554811674544207131412813853349595964490403585112728420124501139008902088625283799686593383080280038118641566198703313983769770057060800*y^290 + 1308111483303858090278066814814637636418637799097773709128168613146826278066783267370809723712720920123258080750760336668929189066102527956484586989983807766133120000*y^289 - 2470269254792234857655193377243597883450123890569836172987525543489294560705128692600863911487350904249605175402880774662609580544202966657981925524528969439375020000*y^288 + 4609160926624535770990787642905737514242304332404694322769407416510513021803471828877221688506886443294995022398058030772917875893451876813063836649426003953955830000*y^287 - 8497601698597991505584170452724165446447858686665444330797753647303859305850554299571739091195042797920590514134913785307722149953832688484517512469705300670681455000*y^286 + 15480543164603779612980004558085567367704309579356837531755065591860223689816167692412796267187600774317089806802972762567611720196596252144426289383406989853395829600*y^285 - 27868223267641601429398339030152157683367275940152246900583150145544524288394737963642270454652444948389611058236455778221852353514864906659847079116591848871927365000*y^284 + 49577327483520458951496306311912190571952449134777752436694341633801983780639273640427527715343797020296070041925310664886803191058567893698115362720104612609294490000*y^283 - 87161504744197076275868259529173410497626314282997659519095296565650002893179039742338579161255655940950849530595452992373777945267714721446169383629312058093218935000*y^282 + 151443519583161138744154804574192415371453337909673311636938723177114935711567044566387613614648909392913591642048040489006800028979523236200514839182135021211312360000*y^281 - 260062177239750599921212500521571531062867926377111225705431996300256814602474297085946729823822099407497728758650273973066677160875392401719884093284432972535648058200*y^280 + 441387418675941878424153622011100743212331184342345639242826812189978359171393230348275600902755197181184211787293844714987718984829624464230458981855299479003155570000*y^279 - 740447502143165534695100537307967168633593753579184089609395784740584659607716349319478979606100202255501771656303918308178342531632346180857749693055635640733881545000*y^278 + 1227764465961450432980543411728357427318362017392952554172603526918778267790914379879293037287528258770541628594756646210578560744649835984081190983745053352489832780000*y^277 - 2012318888937286043566457855960467799821480692619420906417973547313258750402067179971589690637120461454303176530333207097147568581556590471354426229457439520154524955000*y^276 + 3260270208216939071519516239568575051471809861891662297629286181486021501690372530335269572550156516978535452967689977295887707476446074993278017934769798636601523340080*y^275 - 5221546283931850314496818018411388448253306289914236285805562004437738272430807069278221380171375847389524283076456420385749647981184989651837173466736069023132972425000*y^274 + 8266976854011279082251372406181122050250088398332976481390965629272316571446969372233919296017086209807033028893891277346018944293793786433244553441015141421962784800000*y^273 - 12939255550579819499392773054193246845361351594048572480733252687794872494681202396447675154828881925309604169221287717851405988010640475222889118905171760454228516450000*y^272 + 20021627899913177035960773065406032634027774053275650898846779681729335299273004446163807312146043618744750362706494402284217383662221055141936668613414605967305625700000*y^271 - 30628743717846745488141166274846301417230074478281788198795210266116261057695006622865813122138919498601324239408803480132926266281330299195548633970013133230853160378000*y^270 + 46324473617153757753108552139478288541617009319395483500749490497716196639687411350975945904543764786294206340604036429187142135181316386724438276157732344907895241800000*y^269 - 69271526684438649177413390067235390808082029546525213943361083224463368368378137122990759193384549886419382417791765409682295926633834136925662079569706827364935947550000*y^268 + 102417240485161545691395887774272711624946464577036873656921953798614430447265013748193443544748935685903583317569240365808241148336246725775505651806370879693517000700000*y^267 - 149718738872427077655896107931412086562995741311221159176676618083915867328154807275485187605066573767612850427326139616353728715079320761081571932179677105479201733450000*y^266 + 216408549953564029529488472385087115611659957554151520041638535059860882344816669353969231170734718550741231047860394409992951875562200091791265331706051013142465781059200*y^265 - 309298543387625691913571285223267570037485603091706373027602737684014995273829262904936708457571006156644632725061357904609382729032619662598016814016462383950230480530000*y^264 + 437115531999882456732767743895680524235485957381154774128663304687715408569898389960909152516983479183154772124151152596027908911832771674458113816406695488384336132380000*y^263 - 610856034203996790493262391989525701919565131131440136005428071279277789569354335858408666737736304953378914823699631224788889620614679350183831708858121902005740802270000*y^262 + 844140607150987653161112760627611474036431448563584388572560501346872805489381928420148554144420195111968911901467192037360652960327953172777298811683722021321761136520000*y^261 - 1153540133409728065688229301253677208105715044760425703798498247386442488116768044169379678046832867618144911844310167064255298242635749607395003048591420730427064703456400*y^260 + 1558838018121254142821931488180644875818533844270845545673646280251949308265902762391053618982206577862357988978797523059804457084642904874858112227826244230306844193860000*y^259 - 2083184238357813854863328232526205692315534871975142947535183597400836123185246196390414344766845017189238638755031074921006148877914508561224126623326480680164414730110000*y^258 + 2753090650820906336683880369164318137034449671715364013607055424590416926429913113409230049686449564163997988232034369417349011805752947530648956585583220687099028519840000*y^257 - 3598213186247577540419756086850740973174520638895087710295836810455620825205015421586518018211218316724659222090113905725077389913872126957620721339338683110578351469805625*y^256 + 4650866926224319832636672965592173430424792166979689840445128677357539701912051305705554269813402185068422272944711934615629442022008835471497214719396195895084802527246800*y^255 - 5945225915008957775005675142473978238261261078125477732901381660117997454821061329431159814508133203512783643521905814594611167450005181060119257090439578921493570013005000*y^254 + 7516172010249942811440347025385049929091531352758252920705069831178308596081786801488434512969146295642863350234984426121953160898586553013495813824189167554679910793780000*y^253 - 9397776761963842946069840852973774936902618530195518350240814136222776102496944403974206547657570956225577437061022620050574359494364270528292226657817982677008606934895000*y^252 + 11621427176440760405372149907195050028862377165996323634745369808334568912678662373430740698230516761143396993881894263216753405784207531507438235006062739048451222607470000*y^251 - 14213639332814855835428614181039938642572988858222558150219118841219086938510332009744346641610370936156255072444495696328446692732583131444406458043415060551113033497441672*y^250 + 17193641230965858416108547661779574493846456740483086744833694828977460369804920900160094161719614525761183374999389965075297204157089964005910942617959864217246133325400000*y^249 - 20570845002790932257499387801546254967756079243472991870128800764696304297474493807316805722604032725661396499290003191023905055128162179023973660364654462812168320272525000*y^248 + 24342364878953295198868564165392502986714249931836750356439286503509912430672535760422316083544159026433324024029697240506020672522972711557576696995999439510524223989350000*y^247 - 28490766648505998727131030541295853245550793463059711252080128012081801822391276219906419356068101074402730400857772598878685233662753864715450149110702283329547270291475000*y^246 + 32982252214270473820584628297217787757202330314883242202408054075162932933215312706338725466318837008485043193463586161525442670498999768094121113794036525642675898784484000*y^245 - 37765485081410504428588596896260660773184290355515935581920352397682607633121502361838093403158432510794097846987892552652622557355232334404578568965345964255077199320025000*y^244 + 42771246095646656082679438766422495893373805346652845408789879440537465664691344525791175933867058123527357060125822611550496252631047894089748567818538055411482780816650000*y^243 - 47913071959101240233640256979564530523654949050606023905480259778382616017970653201048114463345701976193319609297304606580679399802787827471339943680351986355112285327475000*y^242 + 53088972564212334292360874238052328569921866568592509532699832840720414715004255027745972013757206482134935907053702514431360891697351400191451784784754010922018820657600000*y^241 - 58184249907220427614563414334566350866716034286163852530945381083123843089199068171004162803934713609027457969942594772423570790612399578624112312033976473255986865087856800*y^240 + 63075356498888633967860681191903423343530124654076447209278941075276893421155323045782375044204615361631135494232177672653928015688603575153564126777038017793671573031720000*y^239 - 67634641988588958955567460241168625928182102757187656005562812688395113533850996944904231665046923052277453862459329099640378445536483693995671204734500885480651821331620000*y^238 + 71735751698132659638456133338805504160872801540027603477762823944662944103698596863954746636819312687884663054612435908460368206681942502350123557214253724457077824246480000*y^237 - 75259368892849667231245558503989232527831684770677098724812295104713498241923542145952334442694705023538499713979939092582498308681398772882660873311677653967740557067980000*y^236 + 78098942243984136945218795179458897187747036281173567981663368934196305409632997778693231033014817298185484242195352397352277677491221762043907368673508328571346087306289600*y^235 - 80166017709241402316035961666595257950543098699627446475015700156164338142413811645386226074522339850748782344747673076330223167326128963140644870687084445026190469510340000*y^234 + 81394804001535108666680289711760727876210528170805059081928877657025226219397469089368234444763332184451246819682493325814132167438424011625326735657763237868406755063840000*y^233 - 81745643673955518617829773719139006875676952861196460198661329715891714263446596111305166489783863788694571159422504072908158599194624115123711764690770838031632646249460000*y^232 + 81207134822348433514945519610317274287563138942479750974691228729621017858549055952469164075490162446213052245328811687711135156512182586038575323553150476779777358724760000*y^231 - 79796745336244661661419203770316699300524145765546906416047355215374485962016323817360475185906443761189516508330033971207761974827987452523122746253367205279418828334687200*y^230 + 77559873618831971469755248065246545312525229719999106109413544300054085524427755224649567373257425264432689661568845527726941488422751573800471467692383957800815595930320000*y^229 - 74567426107269475150500454391653089851433570021870798351885040679320691897265646740163486635641127662406426069926084694320645899050966548495523171314816371475120369940420000*y^228 + 70912091600370802781608174173560057140560971172734194402308147223956110224918666400384172406524967781682969247864847249843216167566885987537311659589476699314516006506280000*y^227 - 66703584970063643886365941614446540615503096426213587977954278393766272246266027468769856794534969342609917624431373733470304405971576347175671797786336664834900468825580000*y^226 + 62063199901035122788867641956870774257830436571908109366341739821998765961812339382404440324197027849629288501272101090585142176653226217199053132836241571473052656093713920*y^225 - 57118046374889163234191833443618799934603436799244619631163659440239374760339399382530264532658423045511248007142599580258740959370811104670267322022882663536636320439650000*y^224 + 51995351632746637383591669009213526397867702557159900381731941104791897203358287330285397758428968063851046392152231904899212622207509660305220925428902155596130820848650000*y^223 - 46817174780369988414824791082165164423828084230795682136648381549777496200540525764717368251935563621698863809056643669399867519771822694882022080990561925517974395218225000*y^222 + 41695829682643271857808224627545706916999849985275512588183556493369032861907903060128935703034649468931040156855410782114207710383152441605055071617848032933207405584100000*y^221 - 36730235420437573118378336058265263638702595123392701543590787492795084402898870968422671542036886668540161738175357316244242973964795196286634876725177039920252705464539000*y^220 + 32003326920239615008558143895702191316920069369464333187268088999872708896384721233319539550651839563027724559867634974082906700826655026118724825028460099082505152930150000*y^219 - 27580569381199375752846492425926368062928530886071914394681209937463065919156496939561064553373427887747848779519886322958656450060432384195969927330329640179355428098025000*y^218 + 23509537146046779722856443201767148534578376278254847851575772302510062054256564351219206540118894543613024570617146504027320063003613918884541326131033671785335223168600000*y^217 - 19820447530141197896347353554551453930406332075834669441961064872293461259636228259906775214204555633863586390624296765295134277014423152951258447740148176805992070361225000*y^216 + 16527487596296343994684435094232486793689521928165765290952463768859403416872479081098077490244617367623552409819181042242379875270445501679542393073738442501294224997025200*y^215 - 13630739351806686269852381433205409146002133965724659964117232812579330827487589980051943529232989285266320773197922710595350957251570515926948584049167044397949798478855000*y^214 + 11118496450774124048871361009152071756335492345583430160848287147080467501545680703293427464834309300323218376241209860341524594331360545413754183651502965420686239190580000*y^213 - 8969771314119035524973621685954800151227157699282254664803691826830717334003264437625348734816952890554926545327599338721010140128973575058558683478867261279877440711445000*y^212 + 7156812238920735899769864863452087547480380634331143811170383609810290774997865297040353520782302529347142395519261027995408703944555999949148300004604990739617852053570000*y^211 - 5647481267558102653915960029647639165515819059903909091969084822505583923750348014071030184773416939011655617147150206318978608169582157358246617873552523586888688783085400*y^210 + 4407380744877415551901710804946362148960178390225561980075283365916932685418583482900232596469243488428290721794589281653683310133738107661009113219477047255062478355360000*y^209 - 3401655275293139079460371417669373852008911373225233420658717519806428606213883029548904436391430471106478283536696312212782192471531286131619143473533114022040626584610000*y^208 + 2596432528967710022003568521699328833900521724490757876638054918499592849187408302747472951448531470699630767143951774539225055123197793279206882458059091862233908214340000*y^207 - 1959898284919440888305243958440713703092234299350237846259584400899295489518780657569334404448502444982285044395672613944123079538458797876641995470907753668201323525290000*y^206 + 1463026510665996571116778625980273532231638593354686258476144118636463643465169226751408372156958758930330132095222650910762817283342901941086345329920826292248960111700800*y^205 - 1080003423854940085662011914332385268252493176838703555967834759076146089376120996655190061920067014334430813756288907459564364816890175382818014656124761997867547336430000*y^204 + 788395472716488468518397767001259772920818051016932589327202732390003846636309536920704781700920748264952096964577849166968911727098748717971250776949891321800512759580000*y^203 - 569115766910738050102563943784970187494716865303186661365536876895391843831105329503559936197854727847374777567582936300769010117881281932586554127301786700223170310970000*y^202 + 406242862956881505796929449248470838258190045851793445170582700241691337076074255755737766332332067787831756415005962620979487763505530303113813128857443131158511675920000*y^201 - 286741328688472339126572225910240428185449443759250764055171058954315669106882086894009231301223502079732131138415022593240661235119915121623426564510517592413767159903440*y^200 + 200126555477716596263660124169626206159582247179823257994954675428751862860749641885824421622852807146658382983260066019849707729704016695717076050049216633454611362300000*y^199 - 138108097017869357068700515331973808100306596172064875082818998304192186512221003387715401143676595091597256291873672614498991851544332705970000045617681943829136561550000*y^198 + 94237406860048517662674053492482663637390781312024148917871930477950166105204775375439205400137357205890327350092350912664240143847043394396332169160819440273110832200000*y^197 - 63577976545991636916256125130202344988239670953660127865824555836151310694264865578635354328174860854658885232767784691078271603896806673616703552481785718266448061450000*y^196 + 42409029763172416829733922534402571193087795431467091352177983445391746079885393785506510476807734690371064331489112510953282522543380414270759647035402799858104375908000*y^195 - 27968356499363252457450358251159973421760447259988151513337511823446648794345279406754596812327720725153749218007525685665066551404442814141571325743481227923401157550000*y^194 + 18235634765583994320387742335443332341649158723218713547976397570363758127219103895511749441002362091628043181493213778526333827662664296687180519907806833386270240800000*y^193 - 11754607714297287445076170140038782915199344441379700934477779254121575125026159576145103999024294454569375087958501659976370686417738040084769599423943404209353334637500*y^192 + 7490595192656162231267686581385970862939298028463354671783895231497069563591838413519542712916005142328427505565178320343954992391857000861797829999401481216431593950000*y^191 - 4718812143471885498251931399375569153446040589580037536042387894256083752446591820222100697112068642995248681926919703980890102926435108542902394379272076974028520623800*y^190 + 2938624409241751047747843218227984723548828947952170212546882005600314365666510802498194468866906070695172892710962778618056766817511211481977892271588641838955814220000*y^189 - 1808993388677243457104311303450824507272267440520467327158570090194157341372783421417894528343858125954468250642274420214522666629085130678672818630305211588508955745000*y^188 + 1100772222437433254787368713217191833583631355588305758884556300309532274710936604870669285189817412978767339681564113848521814030982846194002430024970410203871759730000*y^187 - 662081128013912817636230433016651114818496776513518458864240218284967538757089845522413401354132795771776317680940767072481578477244436696345688472782949218504733355000*y^186 + 393607619606109016096970271059675830797854158841033206473086264670658962641582719913865320814433447751888454337416630878154970838081674986386851600646709363837737601600*y^185 - 231280592406892826016888587420369938692273847568577177010759680504180961602990975669320214270197519235916952453345105312015819609474163203110568857664021418262928115000*y^184 + 134314449933680740490682235995775419455358582665602950734343571756903283146866061242292236740295600958186555014360206250543198800904347275671146263297628417389550490000*y^183 - 77090020207891702655348652122368955669574412665615955712002740012999787966598430990335111777155464597493339680576169880016292874008652274266849712361481348178188685000*y^182 + 43726935955857540364359824776592335812576278270736122208649804648257706618254726970576822149841405333384804091376538348149204797780230019289483814746475608174736860000*y^181 - 24510860010622651053885007216968814998865156833981659425278853853778687508875598846320260756380367892668725527438040207398718556977944065895341624857061067267876730200*y^180 + 13577173663681707098307668991503415182592668069197302115129371749323484700967924746801773073082582347515580374539582919541947122272604183835597813219752603306228470000*y^179 - 7431624640119335765586269963561498126896299682023299273229325441227624928413463004987430475289193310256536918136194690066446238097869037333423152350479006375198445000*y^178 + 4019421160339907564998280749413311810359523618554309705200041562287665725893090754062422635640435369626476217604687362162890164953279403375699404243301374003050880000*y^177 - 2147978515197812432331598349326308391619240476927594329193413234087686331598297524262295512232820147247849916687328052825684861591497901019454804773325684343252680000*y^176 + 1134132656024444964271083928444290830774958971817769805814122187598298383083901092810492030458929037746864756010909211891961606920310891738272136920315961333237415040*y^175 - 591623786200570744346591332444290120742171686662701256788856842984517983735401736186138054234966201667207195539148161638209094085510582266769280181521239427781000000*y^174 + 304898361703136539913582730349639412533633409244600447056713687580231214280681886493864105913477064137287117253921099113982692733697423278742033850773951700516000000*y^173 - 155228544326754771576732867238730430168034867339359099265953932386346540546430994896711032389053460804319374295295891186870459342529128862770402390434943971559000000*y^172 + 78068156795911756465491383523572029207315781234999313081123030322957909280661202111796191727828056310944246721610799076437774874020497556714939213902954278094000000*y^171 - 38782910738636498271968360539885446730577931081233392243973988485289846360772171257987395030984161179585772122853518455956158461425164824633584808520192724322760000*y^170 + 19030436789055288515614435825604058144604555589001808865566685069993574487895133825538626512533720307489049677832582299069265228582127349365592199302735641896000000*y^169 - 9223084496450041144628926134777405372714049967476753858092625878220570091721106525096571182609544271831315742984650631785762314729495930284640079048036352761000000*y^168 + 4414680970099131148617052146412673462156734753603901112681770702381798857887753905580317440700492699555341010453151297192412404075659825235723027985519428654000000*y^167 - 2086864360448224791039545698184478317537147151271752645296586774902936272421464739950980366619415659181589200252676324474828130914847760325175212750013647759000000*y^166 + 974173578782812523060787683981583381393825588967002883341359392424992505739338788666571586415114374138198342789467522523328254276274144795674518708224552757510400*y^165 - 449056060086577291445885878523723269408897107925587399460749549880095053843259972669035950156126597251781011428193232919059217396081573644396171229437394047560000*y^164 + 204390692846885019010426584498966068869138677811295696278074145110484397863803399771650017628759929453191590945113951763918860853957082015182364819007730039760000*y^163 - 91852636791752675673475022199933050858846076870303406766843499702093693560078925008393638674282172951699017069683645087093678282072792993647684028724524950040000*y^162 + 40753514988740659935920448055841572719296307984906954868137341873986109843841933177294815378935192107292682684803150986947417683521513084480451386931524935040000*y^161 - 17850638881465302296932213902106488874180007843094906628639275407598180613218117348319501118552716131113419908336203708032776995642474517959562419996403973384800*y^160 + 7718436623275091401975406954527176793080978234042608998232391096052672434960663947779658502518863107006588796843484307718574958807673730249810794809765868995000*y^159 - 3294297027975683620342702505267953036518899516677152552538185738208364920708842704993578818577412611622203379749160876153046904671619795571084836843029411307500*y^158 + 1387791549039082242075936294620361196246445004507453985990016329353368365352680884758071124282944966294407704341077400648683957470906043342344138050892556030000*y^157 - 577009955259566938536621468501322832782497515857853665742161279002258072704674937705445916662439176536475473349810891334948069536524864422038391622535987742500*y^156 + 236760213900054563167284679978607304290108658119738665427106176416410409187208555084299098707942784694966710355148210896146433693696654046720269078666379486600*y^155 - 95866857770398268785081309813699404819662596678577568289449993063071090449366055480842111545797104565621316888840054046245290549955541737830829956095320847500*y^154 + 38302738962498095915328599959492333173469943331938789944714640744852232316445810122765236205096203829840617953691121728400553146936285897148594835524189460000*y^153 - 15099381211477073582391019780909537684628598083050415017479438250403400267031957217446902075147865623235398697523974067387040150154315825824612417776308327500*y^152 + 5872486612049481959196025754646334343329364456141221012216566592270076202590393956936001649937098851661690282119498643908505681516970513836664432199575665000*y^151 - 2253114297367714449300593428531100504908582070958736142322386342105935400395121968232225438067690429853897823713320096056569060347555152365357717118886250740*y^150 + 852720251817018639034891795629169591099276658030505175728813522016173765098824986715531653908648040768241089864757567766274812161611931308481751949924150000*y^149 - 318313504265432785312371922191666176695653152713027174996855153944239781581454787403603791547891164322611709601048660849284334299927563486471665224550337500*y^148 + 117190046346477741856296130558125856096907130849572193779936723342655441975261961531675027734049483382454062241182094044015128050719600487556732470729975000*y^147 - 42547559985090037203466072353931776291841899953762028596733659371159646257708924720684275104644588876875604050950131867739769589937637708205720025974912500*y^146 + 15232352509604647801884558623568525581263245715247218207198441117399070292415026536018155807578507220901365680156246443377793373501045469175258541099984000*y^145 - 5376836198744465546098805735948555757531293812206089350250064739446806920818393566360447320863928274085596357219974752909212374948157578359193993533525000*y^144 + 1871165538518230296996951728758732279256056999229268445212486536228711910741407222188785863204805972799841430772562172341091595724451382678682895322650000*y^143 - 641918962409171560639799335519846787954039674785895915309313791603612436179899258316273621079314725075197110556784407211581724488971752005463649200225000*y^142 + 217063449938225844287372224445599172693264334819147492819818692693213047218879844454206103064144731800548093707973067655746328330587745721369185514100000*y^141 - 72341327952139631377955143165233324281227913767908609879768666128483548282219409979010870348468598798255391957584478729633276332721334164958130372244600*y^140 + 23759272552313610555787582209413509695117081468972174029952173736286432810863428363453679657770418748254038647314264205853906834331373435472642130390000*y^139 - 7689163416121793443877156372957006446876299242731511665477380350850192547041057716849408021929520492195121438443069180416959175626888002151207445065000*y^138 + 2451746062211391576496460300194031567790415848441198279303610433774214627192842643318258967030197365363061354525353991676299748695397991338661534760000*y^137 - 770144458030191234457656291324390213380085560235363754243615137048146324623431674723548614485371136352580863498516082299963543767142332352032121885000*y^136 + 238296549765108906831130899029895977134748166997693504487654038172357656423165525597000650661928592877772110568218097846041629838637479343634066073200*y^135 - 72620515975849171559123255600150218771478069163927682689643333723583260940559420147277379995369166709312024453452559430162205938236709923731898555000*y^134 + 21794286989392584521662970325688979429738099144381432206970227814020508322586701679138116280805080266624437349972324666610830001620871511265564180000*y^133 - 6440387998605437962700437129919708853326878585898560609614100675633697167069134806617562523391350280096026445720754149125092418564428756099704545000*y^132 + 1873743675328063125155078257559695963563330120787173960823111109371549783258690332943460293333726895811186707033691079188528124356020235889237770000*y^131 - 536634119808531115277207432954546123329837946333946097119946882092888387318703458106965105872177492346896144354507315553831900878319317760342833400*y^130 + 151271069713468954271235358126721951608129089875671908983776429060715542584553476562921805742685691993487285230305092474653108069998398241111440000*y^129 - 41964380211508115009205126705940314918293868546465930543708231395844732536998706436540073016215970358929358770199237734352946002323434137499093125*y^128 + 11454843942248934280727908602146437667985727897198049282272063163170215705637442176903327017549766187214313155120001901240699171237892782939385000*y^127 - 3076204726513218221871758615378426263394091094268680285915815501528026017065649119787379439353522369000076024164291168018696028877227262161922500*y^126 + 812627212719740074279833393127277514434808726033569170149788944624341438328845681960936427896533247049509738135096613101987122000781469309312272*y^125 - 211129035562785810785207794606147508530826177989267015024783043893504073393551941813096635947511339959238271672113145655446437999039083105387500*y^124 + 53940306000783010350959648159761931681670035895687552340667368603597098158459533858920111689897651616715958928146514771532924464839698924575000*y^123 - 13549208515107442789158694142743310708043284627857061838159578769420308408348015288962224664678314045389486801763205844950319721641460768862500*y^122 + 3345596435506374809129419501598280802255923207148378905560064651447348880050850798286506478173225124380200485373812338844201057932720284300000*y^121 - 811921440040544236870927766616690123507927578322365061516003260468479497980702237516366586270795961682211931351593836243235347138418078034500*y^120 + 193621965033516113721206939574727374445451727739832685576153400747650118755334396863998391002574553660304912087025557132727030319177602075000*y^119 - 45364091619306078005045659443248486890703764665632773533806870214545315445416045558702584846263187242617235437612476448576774694708876637500*y^118 + 10440012366010209657330533223559509516750043242042028762555679647837767771869969298543840440727809741144974965622741668771296969837206050000*y^117 - 2359568012257929693916908370991898140404901077567969769048378870931849312757717948658896471423864045693560683759562678511653700454211112500*y^116 + 523618918894107528597908709806202091679418047822039900053170685097224299665886626519782938702060960400866683909078604831455673352969283400*y^115 - 114065728820115806267685714454781157408161389155147993614776792886336528530701828211091026018473674855641384339347280151399753987102797500*y^114 + 24386788884474152333383196760888569549608670669752088185354611759324600368047838557935888349961435027450024618762168316961044514150560000*y^113 - 5115792965245505965975806746250758092399343661538340281457434026244455151465480984616501330839682224198984124851717239758535946965990000*y^112 + 1052752084834779321590090150722062120484265125228847644453544982166948096603138334435069500988346176056120962629514121410718160731740000*y^111 - 212464511666655463084545466782088900679551688909821979153351805491874979496269736585986753835829864622235321548865577230163119711314800*y^110 + 42041907327454759739701387562564812151017603489980963674885228645630954748317451024390255301207564505209363861081513753045409548920000*y^109 - 8154527265739581369466148613067146427289003741225346163384814375076989684715572672754666782545115434099451652509808630288184003170000*y^108 + 1549944131344603037332684637463666865533994399655163501864379088647962840910861367342372640061372832843760321377638101108301510980000*y^107 - 288609307476055868739453001278955370730043104460427495813361667633386009852700234457505365199191332585763393786475385560833208630000*y^106 + 52632407040794704664743257007431000941736893028482261602092622398518566958083827703218182728798763447898356544286263861416464713600*y^105 - 9397481779018965478570749110125959387785387003446854463025132771192330082331625803428302166052833876812175644424102578520902910000*y^104 + 1642278369148945423439548388177352126020747243320809517810217377489921761960866645259314941640301065850671471841105304984235460000*y^103 - 280811802831049710786083275420039004446133533837258007822025204025480263689293084522664920955158914943269296647891807927987890000*y^102 + 46964299748246442139414318371315093636321254445128853928486525892386558573753248417361882472919306453273013477341809587250040000*y^101 - 7679781206451346776511846156528382454625104179265356399543939519735973435346126431677176400571661826692358489580560672505554160*y^100 + 1227429549682161292755377534287236679233011152548484273037965017219021454313087589770677726724789321488997329238678025908700000*y^99 - 191664866590747583569057375227817714896006473347476250837158575243979690815450585784380275939675304144180970969446568714450000*y^98 + 29228798951010718832550278579772240513220983292048508200025738765080738688837436151518019127933366747527306212808035376800000*y^97 - 4351296186574512381060498956543287764491067882329191120817572591105116499556226208685039485509672773053062598752345653175000*y^96 + 632083024997139693248788269476814433157649860801503552287184229023690607303957070314247841052984044927708040660867052777000*y^95 - 89552426985576065110312515808755276123449100381998663950013673330558671247577930522748376764814303442233441831550561762500*y^94 + 12368558004820096349971701241735518000620633112266963003314517895384714628998274248383332875818872605258716203008519450000*y^93 - 1664485915374732982350815495547864331823997570869416645923231830113187299817071568154696061731042924634190540309296987500*y^92 + 218137755291275547544551318504848316991179168831889361004184228447882096272322769615715153389253488698497905567600175000*y^91 - 27824682563820480953460545960396207545097076202112109603644832695352071835625171057648999565651445002875066176845000100*y^90 + 3452431833933746861255185960482055116921123400000123890987847166997568711940808430672654526654208382517068940395740000*y^89 - 416436992479289524444563193940548141188254772080915694523867720206024499553327288812370859814316002770251122189102500*y^88 + 48800297969257925768945168227403735076167786579709212310338297461765696891867341385918424278773223761694692760035000*y^87 - 5551995446491880365103361013784287241471128331097899004275782127992640782318547109679487091021807001091935251097500*y^86 + 612797786179638560084135675532126677347401544138292168172471353592342811481576536597778682132032601617847826645200*y^85 - 65568746003611880170130327820028660088782894324968321116138047652755713715929010084102703903855452086501844492500*y^84 + 6795773645896194286757592402888730338712583328208358437642991148797732558161317056183674058259949599326798805000*y^83 - 681667912404221641816049424723951966517753084888401587955180281130776733201635077158393146032886616588152857500*y^82 + 66115790720362149335599678068886299997877729531572623892747581048538410139760865027856273201385912690445520000*y^81 - 6194654869359304364988462374290055943831006673648464872936163284062684621676848212124891865995523386780175400*y^80 + 560095376976428966092989364764019524758680531071289771513215486804944360006948301277114996925454194103090000*y^79 - 48816043089591370361372541834983212783311918766047510710521824555661855313700013628009984705022276447165000*y^78 + 4096451168357317792562730783355234639159042134213777122561271990684910935415385758993844870351519701860000*y^77 - 330558996235633585263853334055097532014114601654651554021304318763931448328232688749736201112436324335000*y^76 + 25615081433788704489465654435014616598427076583121234146905777799118766348885795802724656211693104975920*y^75 - 1903316069127122965750499408165665614728175025495455101270045707115308191494050373706456780412323525000*y^74 + 135395711718728285955815386882248743561533239016259242964655979165136078991165323387914253064659400000*y^73 - 9205110138311661488492631785523794964359179960653517530679356536210863161085462340796405527031225000*y^72 + 596990438056242703806597148225912677020346423913985087151822009226908256107082064683160266738350000*y^71 - 36857670523472375626320345673078087016038779215559079293721184917487379290089414428264677337759000*y^70 + 2161378015115243607139502552664376385291710266473264262847015953944974910200434790968021624900000*y^69 - 120086734099681487072380396954695602691872866559859698618089418359970423261263350184737936525000*y^68 + 6304090004203557221504271996592357938533727426199531270738610742520320727045219618550318850000*y^67 - 311728684326988437756940427373391321176952137641801950046508536760034961713315322806668975000*y^66 + 14469691896672518253684795222254999347818086037351334472488484167850414047000482676212854400*y^65 - 628024821904189160316180348188151013360159984260040558701757125340729776345507060599516250*y^64 + 25374740278957137790552741340935394479198383202425881159666954559221405104868972145435000*y^63 - 949552786642871966118189489090107286779835321337238265090342492186892569719806642927500*y^62 + 32715905958242879632857876256508281067157427020232715996444142548256370837803840990000*y^61 - 1030627835585961137729935913996105943477118945097471945240094819243193654092318183300*y^60 + 29440901797351775025612338049784020476398370170371508624264753930009340718691245000*y^59 - 755004291238954531952098498234101944990221895125435674658774477776246494814307500*y^58 + 17167597330925948388992144623180598964787952530540142754173838953649827187780000*y^57 - 340721813322694658138075293705185982239872316936207568589104742464498659367500*y^56 + 5781945923051788138100671650754671213767530226796249648784807750912704522600*y^55 - 81565360051258938031070314145764572883249759534675728604484566097311672500*y^54 + 918279879655578578896653177081473352589489630954339486008999597467910000*y^53 - 7735033951326392988163162489719975792632382460359185117367606272377500*y^52 + 43333523536842537748813235236526475028752842915177507660322724215000*y^51 - 121096421938573667133669862852758916518706574721865911817888160820*y^50)*(y - 1)^349*y^58

Worse, the expression has lots of large (more than 100 digits) numbers, so SageMath encounters errors and thinks that the integral is sometimes negative, which is ridiculous.

Where are we integrating? The placebo proportion is represented by $y$, the hydroxychloroquine proportion by $x$. We want to know the chance that the hydroxychloroquine proportion is at least $10\%$ better than the placebo proportion. That is, we need $x \leq 0.9y$.

(Of course, we are assuming that our prior knowledge is uniform.)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.6683714444990453

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.8470851993226352

So there's an 85% chance that hydroxychloroquine is better than the placebo, and a 67% chance that it reduces cases by at least 10%.

Before I figured out how to use SageMath's built-in numerical integration, I did it using Monte Carlo integration, which was much slower and less accurate.

#This takes a long time---1 minute? 2 minutes?
import random
for j in range(10): #Do it ten times to make sure the results don't jump around
    total_volume = 0
    effective_volume = 0
    for i in range(100000):
        x = random.random()
        y = random.random()
        value = f_2d(x,y)
        total_volume += value
        if x <= 0.9*y:
            effective_volume += value
    print(effective_volume/total_volume)

0.6894516820972163
0.6249901715346099
0.6759584915783055
0.6472110465519125
0.6981347683159699
0.687890644602133
0.6761077587198506
0.6821597024531272
0.6786683281125644
0.639106862842741

When I analyze Horby et al. later on, I have to use my own integration method because SageMath won't cooperate. Let's check that my method gives the right answers for this problem.

part_volume =  0
total_volume = 0
density = 100
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        value = f_2d(xx, yy)
        if 0.9 * yy >= xx:
            part_volume+= value
        total_volume += value

print(n(part_volume), n(total_volume), n(part_volume/total_volume))

0.0384530013381997 0.0590597684857099 0.651086218658371

part_volume =  0
total_volume = 0
density = 300
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        value = f_2d(xx, yy)
        if 0.9 * yy >= xx:
            part_volume+= value
        total_volume += value

print(n(part_volume), n(total_volume), n(part_volume/total_volume))

0.356537055127285 0.531537916371368 0.670765046379466

Let's look at the participants who took the drug completely: 43/312 in the experimental group got infected, versus 50/336 in the control group.

(See the paper's supplemental appendix.)

n(43/312 / (50/336))

0.926153846153846

var('x y')
f_2d(x,y) = x^43 * (1-x)^(312-43) * binomial(312,43,hold=True) * y^50 * (1-y)^(336-50) * binomial(336,50,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.43644423926443965

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.6527656841888028

What about the participants who took some of the drug? 4/37 in the experimental group got infected, versus 3/15 in the control group.

n( 4/37/(3/15))

0.540540540540541

var('x y')
f_2d(x,y) = x^4 * (1-x)^(37-4) * binomial(37,4,hold=True) * y^3 * (1-y)^(15-3) * binomial(15,3,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.7860010827253943

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.8296253815322459

What about the participants who took none of the drug? 2/65 in the experimental group got infected, versus 5/56 in the control group.

n( 2/65/(5/56))

0.344615384615385

var('x y')
f_2d(x,y) = x^2 * (1-x)^(65-2) * binomial(65,2,hold=True) * y^5 * (1-y)^(56-5) * binomial(56,5,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.8764540990473931

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.9058717565052852

Let's compare the percentage sick in the various subgroups of this study.

print(n(43/312))
print(n(50/336))

0.137820512820513
0.148809523809524

print(n(4/37))
print(n(3/15))

0.108108108108108
0.200000000000000

print(n(2/65))
print(n(5/56))

0.0307692307692308
0.0892857142857143

Skipper et al.

At 14 days, 24% (49 of 201) of participants receiving hydroxychloroquine had ongoing symptoms compared with 30% (59 of 194) receiving placebo (P = 0.21).

Same procedure as last time.

n( (49/201) / (59/194))

0.801585293869635

ineffective_proportion = 59/194
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5

f(x) = x^49 * (1-x)^(201-49) * binomial(201,49, hold=True)

show(f)

print(n(f(very_effective_proportion)))
print(n(f(effective_proportion)))
print(n(f(ineffective_proportion)))
print(n(f(dangerous_proportion)))
print(n(f(very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

0.000213115038855560
0.0411358667778100
0.0107795231338404
0.00129091347725268
2.42500308651183e-10

total = f(very_effective_proportion) + f(effective_proportion) + f(ineffective_proportion) + f(dangerous_proportion) + f(very_dangerous_proportion)

print(n(f(very_effective_proportion)/total))
print(n(f(effective_proportion)/total))
print(n(f(ineffective_proportion)/total))
print(n(f(dangerous_proportion)/total))
print(n(f(very_dangerous_proportion)/total))

0.00398946757865432
0.770054556971662
0.201790348943686
0.0241656219664439
4.53955349360988e-9

So hydroxychloroquine is probably effective.

var('x y')
f_2d(x,y) = x^49 * (1-x)^(201-49) * binomial(201,49, hold=True) * y^59 * (1-y)^(194-59) * binomial(194,59,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.7559254454549569

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.9100585108660569

So there's a 91% chance that hydroxychloroquine is better than the placebo, and a 76% chance that it reduces cases by at least 10%.

The following is speculative. I might be comparing apples to oranges, and I'm not sure that a uniform prior is a reasonable choice.

Since Skipper et al. and Boulware et al. are kind of measuring the same thing, we can try to combine their statistics. That is, assume that the chance of developing symptoms in Boulware et al. is some fraction $z$ of the chance in Skipper. Then we can "import" Boulware's data into Skipper.

var('x y k')
f_3d(x,y,k) = x^49 * (1-x)^(201-49) * binomial(201,49, hold=True) * y^59 * (1-y)^(194-59) * binomial(194,59,hold=True)  * (k*x)^49 * (1-k*x)^(414-49) * binomial(414,49,hold=True) * (k*y)^58 * (1-k*y)^(407-58) * binomial(407,58,hold=True)
show(f_3d)

from scipy.integrate import tplquad

#The ordering of the variables is weird
gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: 0.9*y 
effective_volume = tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(effective_volume)

5.411877197980955e-09

gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: y 
safe_volume = tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(safe_volume)

6.528422087537368e-09

gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: 1
total_volume= tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(total_volume)

6.843994100902817e-09

print(safe_volume/total_volume)
print(effective_volume/total_volume)

0.9538906654925636
0.7907483726888447

So it looks like a 95% chance that hydroxychloroquine does not make patients more sick, and a 79% chance that it reduces cases by at least 10%. I don't trust these numbers. Even if it was valid to compare the studies, I'm not sure that this statistical method is the right one. It gives probabilites that are lower than I expected.

Let's do the integrals a slower way.

effective_volume =  0
safe_volume =  0
total_volume = 0
density = 30
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        for k in range(density):
            kk = k / density
            value = f_3d(xx, yy,kk)
            if  yy >= xx:
                safe_volume+= value
            if  0.9 * yy >= xx:
                effective_volume+= value
            total_volume += value
effective_volume = effective_volume / density ^ 3
safe_volume = safe_volume / density ^ 3
total_volume = total_volume / density ^ 3

print(n(effective_volume))
print(n(safe_volume))
print(n(total_volume))

print(n(safe_volume/total_volume))
print(n(effective_volume/total_volume))

6.09161669138185e-9
6.76651578836300e-9
6.84282719741760e-9
0.988847970750541
0.890219278616410

effective_volume =  0
safe_volume =  0
total_volume = 0
density = 40
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        for k in range(density):
            kk = k / density
            value = f_3d(xx, yy,kk)
            if  yy >= xx:
                safe_volume+= value
            if  0.9 * yy >= xx:
                effective_volume+= value
            total_volume += value
effective_volume = effective_volume / density ^ 3
safe_volume = safe_volume / density ^ 3
total_volume = total_volume / density ^ 3

print(n(effective_volume))
print(n(safe_volume))
print(n(total_volume))

print(n(safe_volume/total_volume))
print(n(effective_volume/total_volume))

5.06957969207133e-9
6.72430911409873e-9
6.84272611467702e-9
0.982694470216440
0.740871343834375

effective_volume =  0
safe_volume =  0
total_volume = 0
density = 50
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        for k in range(density):
            kk = k / density
            value = f_3d(xx, yy,kk)
            if  yy >= xx:
                safe_volume+= value
            if  0.9 * yy >= xx:
                effective_volume+= value
            total_volume += value
effective_volume = effective_volume / density ^ 3
safe_volume = safe_volume / density ^ 3
total_volume = total_volume / density ^ 3

print(n(effective_volume))
print(n(safe_volume))
print(n(total_volume))

print(n(safe_volume/total_volume))
print(n(effective_volume/total_volume))

5.36078489855872e-9
6.69338517816548e-9
6.84272621368508e-9
0.978175213963562
0.783428231840906

What is the constant $k$?

from scipy.integrate import dblquad

density = 100
points = [ (i/density, dblquad( f_3d, 0, 1, gfun, hfun, args = [i/density] )[0] ) for i in range(density+1)]
scatter_plot(points)

The constant $k$ is somewhere around 0.5, which makes sense because all the proportions from Boulware are half of their counterparts in Skipper, as you can see:

n(49/414 / (49/201))

0.485507246376812

n( 58/407 / (59/194))

0.468579519426977

Skipper et al. also looks at hospitalization/death rates:

The incidence of hospitalization or death was 3.2% (15 of 465) among participants with known vital status. With hydroxychloroquine, 4 hospitalizations and 1 nonhospitalized death occurred (n = 5 events). With placebo, 10 hospitalizations and 1 hospitalized death occurred (n = 10 events); of these hospitalizations, 2 were not COVID-19–related (nonstudy medicine overdose and syncope). The incidence of hospitalization or death did not differ between groups (P = 0.29).

Based on the flow chart (Figure 1), the placebo group had 234 patients, and the hydroxychloroquine group had 231 patients. This does add up to 465. I'm going to decrease the placebo group's hospitalizations by 2 because those 2 patients were not hospitalized for COVID-19.

n(8/234)

0.0341880341880342

n(5/235)

0.0212765957446809

n(5/235 / (8/234))

0.622340425531915

n(5/235 / (8/234)) ^ -1

1.60683760683761

var('x y')
f_2d(x,y) = x^5 * (1-x)^(235-5) * binomial(235,5, hold=True) * y^8 * (1-y)^(234-8) * binomial(234,8,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.7314235205271544

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.7939248539228614

So there's a 73% chance that hydroxychloroquine decreases hospitalizations/deaths by at least 10%, and a 79% chance that it decreases hospitalizations. I'm surprised that we can get such high probabilities from just a few hospitalizations, but the placebo group did have 60% more of them.

Mitjà et al.

The clinical outcome of risk of hospitalization was similar in the control arm (7.1%, 11/157) and the intervention arm (5.9%, 8/136;RR 0.75 [95% CI 0.32; 1.77]) (Table 2).

Here we go again.

n( (8/136) / (11/157))

0.839572192513369

ineffective_proportion = 11/157
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5

f(x) = x^8 * (1-x)^(136-8) * binomial(136,8, hold=True)

show(f)

print(n(f(very_effective_proportion)))
print(n(f(effective_proportion)))
print(n(f(ineffective_proportion)))
print(n(f(dangerous_proportion)))
print(n(f(very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

0.0555971684750819
0.140876262292071
0.125210622945351
0.101946062708870
0.0235383063646503

total = f(very_effective_proportion) + f(effective_proportion) + f(ineffective_proportion) + f(dangerous_proportion) + f(very_dangerous_proportion)

print(n(f(very_effective_proportion)/total))
print(n(f(effective_proportion)/total))
print(n(f(ineffective_proportion)/total))
print(n(f(dangerous_proportion)/total))
print(n(f(very_dangerous_proportion)/total))

0.124331606710266
0.315040720930965
0.280007747786041
0.227981354483190
0.0526385700895380

This is what inconclusive data looks like. The results lean toward hydroxychloroquine reducing hospitalizations, but there's still a 56% chance that it's ineffective or harmful.

var('x y')
f_2d(x,y) = x^8 * (1-x)^(136-8) * binomial(136,8, hold=True) * y^11 * (1-y)^(157-11) * binomial(157,11,hold=True)
show(f_2d)

def inner1(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=0.9y'''
    return numerical_integral(f_2d, 0, 0.9*y, params = [y])[0]

def inner2(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=1'''
    return numerical_integral(f_2d, 0, 1, params = [y])[0]

def inner3(y):
    '''For constant y, this function does the integral with
    respect to x, from x=0 to x=y'''
    return numerical_integral(f_2d, 0, y, params = [y])[0]

part_volume =  numerical_integral(inner1, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.5456222630855441

part_volume =  numerical_integral(inner3, 0, 1)[0]
total_volume = numerical_integral(inner2, 0, 1)[0]
print(part_volume/total_volume)

0.6407564605821422

So there's a 64% chance that hydroxychloroquine is better than the placebo, and only a 54% chance that it reduces cases by at least 10%. The prior probabilities were 50% and 45%, so we're not updating by very much.

Let's bring in the hospitalization data from Skipper et al. and use my highly speculative method.

var('x y k')
f_3d(x,y,k) = x^8 * (1-x)^(136-8) * binomial(136,8, hold=True) * y^11 * (1-y)^(157-11) * binomial(157,11,hold=True) * (k*x)^5 * (1-(k*x))^(235-5) * binomial(235,5, hold=True) * (k*y)^8 * (1-(k*y))^(234-8) * binomial(234,8,hold=True)
show(f_3d)

from scipy.integrate import tplquad

#The ordering of the variables is weird
gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: 0.9*y 
effective_volume = tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(effective_volume)

1.624275021015414e-07

gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: y 
safe_volume = tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(safe_volume)

1.8566276539272728e-07

gfun = lambda k: 0
hfun = lambda k: 1
qfun = lambda k,y: 0
rfun = lambda k,y: 1
total_volume= tplquad( f_3d, 0, 1, gfun, hfun, qfun, rfun )[0]
print(total_volume)

2.3381682205960883e-07

print(safe_volume/total_volume)
print(effective_volume/total_volume)

0.7940522147093196
0.6946784267734699

So we have a 79% chance that hydroxychloroquine is safer than the placebo, and a 69% chance that it reduces cases by at least 10%. If possible, I trust these numbers even less than when I combined two studies up above. The two following expressions show that it is not the case that the results of Mitjà et al. are just a multiple of Skipper et al. That breaks the critical assumptions of my methods. It is much better to do the statistics separately for each study and then make an intuitive judgement about what they mean together.

(I'm sorry I'm using so much bold, but I want to be very clear about which sections of this notebook are unreliable.)

n(  (5/235) /(8/136)  )

0.361702127659574

n( 8/234 / (11/157))

0.487956487956488

And thus the constant $k$ has a wide distribution, indicating that there really isn't a value of $k$ that works:

from scipy.integrate import dblquad

density = 100
points = [ (i/density, dblquad( f_3d, 0, 1, gfun, hfun, args = [i/density] )[0] ) for i in range(density+1)]
scatter_plot(points)

Horby et al

Let's do the big British randomized control trial. This one studied whether hydroxychloroquine prevented deaths, and the answer was no.

Results: 1561 patients randomly allocated to receive hydroxychloroquine were compared with 3155 patients concurrently allocated to usual care. Overall, 418 (26.8%) patients allocated hydroxychloroquine and 788 (25.0%) patients allocated usual care died within 28 days (rate ratio 1.09; 95% confidence interval [CI] 0.96 to 1.23; P=0.18).

n((418/1561) / (788/3155))

1.07212771976834

ineffective_proportion = 788/3155
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5

f(x) = x^418 * (1-x)^(1561-418) * binomial(1561,418,hold=True) 

show(f)

print(n(f(very_effective_proportion)))
print(n(f(effective_proportion)))
print(n(f(ineffective_proportion)))
print(n(f(dangerous_proportion)))
print(n(f(very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

2.39677198943557e-52
8.12393127425579e-6
0.00602380538451211
0.0188375797027912
1.01781031504652e-19

total = f(very_effective_proportion) + f(effective_proportion) + f(ineffective_proportion) + f(dangerous_proportion) + f(very_dangerous_proportion)

print(n(f(very_effective_proportion)/total))
print(n(f(effective_proportion)/total))
print(n(f(ineffective_proportion)/total))
print(n(f(dangerous_proportion)/total))
print(n(f(very_dangerous_proportion)/total))

9.63739166561420e-51
0.000326662310389288
0.242216498122754
0.757456839566857
4.09260317236747e-18

So hydroxychloroquine is probably dangerous.

var('x y')
f_2d(x,y) = x^418 * (1-x)^(1561-418) * y^788 * (1-y)^(3155-788)  
show(f_2d)
#The binomials are just constants, so I could leave them out.

dangerous_volume =  0
unsafe_volume =  0
safe_volume =  0
effective_volume =  0
total_volume = 0
density = 100
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        value = f_2d(xx, yy)
        if 1.1 * yy <= xx: #Raises death rates by at least 10%
            dangerous_volume+= value
        if  yy < xx: #Raises death rates
            unsafe_volume+= value
        if xx <= yy: #Does not raise death rates
            safe_volume+= value
        if  0.9 * yy >= xx: #Lowers death rates by at least 10% 
            effective_volume+= value
        total_volume += value

print(n(dangerous_volume/total_volume))
print(n(unsafe_volume/total_volume))
print(n(safe_volume/total_volume))
print(n(effective_volume/total_volume))

0.301958885090088
0.839044451427119
0.160955548572881
0.000510455405632987

dangerous_volume =  0
unsafe_volume =  0
safe_volume =  0
effective_volume =  0
total_volume = 0
density = 300
for x in range(density):
    xx = x / density
    for y in range(density):
        yy = y / density
        value = f_2d(xx, yy)
        if 1.1 * yy <= xx: #Raises death rates by at least 10%
            dangerous_volume+= value
        if  yy < xx: #Raises death rates
            unsafe_volume+= value
        if xx <= yy: #Does not raise death rates
            safe_volume+= value
        if  0.9 * yy >= xx: #Lowers death rates by at least 10% 
            effective_volume+= value
        total_volume += value

print(n(dangerous_volume/total_volume))
print(n(unsafe_volume/total_volume))
print(n(safe_volume/total_volume))
print(n(effective_volume/total_volume))

0.306166791525411
0.888306834715452
0.111693165284547
0.000524937366943429

Ah, now the result is more subtle. Assume we start with a uniform prior. The data tells us that there's an $89\%$ chance that hydroxychloroquine increases the death rate, but there's only a $31\%$ chance that it raises the death rate by at least $10\%$.

Of course, the chance that hydroxychloroquine lowers the death rate by at least $10\%$ is miniscule: $0.05\%$.