Vitamin_D_statistics.ipynb

Kernel: SageMath 9.1

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

Castillo et al

Of 50 patients treated with calcifediol, one required admission to the ICU (2%), while of 26 untreated patients, 13 required admission (50 %)

n((1/50) / (13/26))

0.0400000000000000

ineffective_proportion = 13/26
effective_proportion = ineffective_proportion * 0.9
very_effective_proportion = ineffective_proportion * 0.5
very_very_effective_proportion = ineffective_proportion * 0.1
dangerous_proportion = ineffective_proportion * 1.1
very_dangerous_proportion = ineffective_proportion * 1.5
very_very_dangerous_proportion = ineffective_proportion * 1.9

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

show(f)

print(n(f(very_very_effective_proportion)))
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)))
print(n(f(very_very_dangerous_proportion)))
#The n() function makes the output a decimal instead of a huge fraction

0.202486777043982
9.43869427378229e-6
4.26534636264292e-12
4.44089209850063e-14
2.79734395057776e-16
1.18329135783152e-28
8.43769498715155e-63

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

print(n(f(very_very_effective_proportion)/total))
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))
print(n(f(very_very_dangerous_proportion)/total))

0.999953388271731
0.0000466117069850912
2.10638324622593e-11
2.19307411855358e-13
1.38143023577983e-15
5.84352331470649e-28
4.16684082525198e-62

Well, that's decisive

var('x y')
f_2d(x,y) = x^1 * (1-x)^(50-1) * y^13 * (1-y)^(26-13)  
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
very_very_effective_volume = 0
very_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
        if  0.1 * yy >= xx: #Lowers death rates by at least 90% 
            very_very_effective_volume+= value
        if  0.5 * yy >= xx: #Lowers death rates by at least 50% 
            very_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))
print(n(very_effective_volume/total_volume))
print(n(very_very_effective_volume/total_volume))

8.30296275867362e-8
1.80527874248038e-7
0.999999819472126
0.999999267348029
0.999830190538604
0.715973398225514

dangerous_volume =  0
unsafe_volume =  0
safe_volume =  0
effective_volume =  0
very_very_effective_volume = 0
very_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
        if  0.1 * yy >= xx: #Lowers death rates by at least 90% 
            very_very_effective_volume+= value
        if  0.5 * yy >= xx: #Lowers death rates by at least 50% 
            very_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))
print(n(very_effective_volume/total_volume))
print(n(very_very_effective_volume/total_volume))

8.27504631854543e-8
2.16833411622758e-7
0.999999783166588
0.999999265181577
0.999817480443293
0.717983321315720

Assuming the study's data is correct, Calcifediol certainly reduces ICU rates by at least 50% (a relative reduction). And it probably reduces ICU rates by at least 90%.

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

The graph shows that the control group's ICU rate is somewhere around 50%, but it could be anywhere in the bell curve, like 60% or 40%. But the experimental group's ICU rate is only a few percent.