Attribute VB_Name = "test_wilcoxon_os"
'Created by Peter Stikker
'Companion website: https://PeterStatistics.com
'YouTube channel: https://www.youtube.com/stikpet
'Donations welcome at Patreon: https://www.patreon.com/bePatron?u=19398076

Public Sub ts_wilcoxon_os_addHelp()
Application.MacroOptions _
    Macro:="ts_wilcoxon_os", _
    Description:="one-sample Wilcoxon signed rank test", _
    category:=14, _
    ArgumentDescriptions:=Array( _
        "range with the data as numbers", _
        "optional hypothesized median, otherwise the midrange will be used", _
        "optional boolean to use a tie correction (default is True)", _
        "optional method to use for approximation. Either " & Chr(34) & "wilcoxon" & Chr(34) & " (default), " & Chr(34) & "exact" & Chr(34) & ", " & Chr(34) & "imanz" & Chr(34) & " or " & Chr(34) & "imant" & Chr(34) & " for Iman's z or t approximation", _
        "optional method to deal with scores equal to hypMed. Either " & Chr(34) & "wilcoxon" & Chr(34) & " (default), " & Chr(34) & "pratt" & Chr(34) & " or " & Chr(34) & "zsplit" & Chr(34), _
        "optional boolean to use a continuity correction (default is False)", _
        "output to show, either " & Chr(34) & "pvalue (default)" & Chr(34) & ", " & Chr(34) & "statistic" & Chr(34) & ", or " & Chr(34) & "w " & Chr(34))
        
Application.MacroOptions _
    Macro:="ts_wilcoxon_os_arr", _
    Description:="one-sample Wilcoxon signed rank test" & vbNewLine & "array function, requires 2 rows and 5 columns as output", _
    category:=14, _
    ArgumentDescriptions:=Array( _
        "range with the data as numbers", _
        "optional hypothesized median, otherwise the midrange will be used", _
        "optional boolean to use a tie correction (default is True)", _
        "optional method to use for approximation. Either " & Chr(34) & "wilcoxon" & Chr(34) & " (default), " & Chr(34) & "imanz" & Chr(34) & " or " & Chr(34) & "imant" & Chr(34) & " for Iman's z or t approximation", _
        "optional method to deal with scores equal to hypMed. Either " & Chr(34) & "wilcoxon" & Chr(34) & " (default), " & Chr(34) & "pratt" & Chr(34) & " or " & Chr(34) & "zsplit" & Chr(34), _
        "optional boolean to use a continuity correction (default is False)")
        
End Sub

Function ts_wilcoxon_os(data As Range, _
                    Optional hypMed = "none", _
                    Optional ties = True, _
                    Optional appr = "wilcoxon", _
                    Optional eqMed = "wilcoxon", _
                    Optional cc = False, _
                    Optional out = "pvalue")
Attribute ts_wilcoxon_os.VB_Description = "perform a one-sample Wilcoxon Signed Rank Test"
Attribute ts_wilcoxon_os.VB_ProcData.VB_Invoke_Func = " \n14"

If hypMed = "none" Then
    hypMed = (WorksheetFunction.Min(data) + WorksheetFunction.Max(data)) / 2
End If

Dim n, nr As Integer
n = WorksheetFunction.Count(data)

If eqMed = "wilcoxon" Or appr = "exact" Then
    nr = n - WorksheetFunction.CountIf(data, hypMed)
Else
    nr = n
End If

Dim absDiffs() As Double
ReDim absDiffs(0 To nr - 1)
Dim scores() As Double
ReDim scores(0 To nr - 1)
Dim i, k As Integer
k = 0
i = 1
Do While k < nr
    If WorksheetFunction.IsNumber(data(i)) Then
        
        If eqMed = "wilcoxon" Then
            If data(i) <> hypMed Then
                absDiffs(k) = Abs(data(i) - hypMed)
                scores(k) = data(i)
                k = k + 1
            End If
        
        Else
            absDiffs(k) = Abs(data(i) - hypMed)
            scores(k) = data(i)
            k = k + 1
        End If
    End If
    
    i = i + 1
Loop
    

'sort scores based on absolute differences
changes = 1
Do While changes <> 0
    changes = 0
    For i = 1 To nr - 1
        If absDiffs(i - 1) > absDiffs(i) Then
            ff1 = absDiffs(i)
            ff2 = scores(i)
            absDiffs(i) = absDiffs(i - 1)
            scores(i) = scores(i - 1)
            absDiffs(i - 1) = ff1
            scores(i - 1) = ff2
            
            changes = 1
        End If
    Next i
Loop

'we need the ranks for which we need the rank frequencies
'store for each score how often it occurs
'also check if ties actually occur
maxRankFreq = 1
Dim Rfreq As Variant
ReDim Rfreq(1 To nr, 1 To 3)
For i = 0 To nr - 1
    For j = 0 To nr - 1
        If absDiffs(j) = absDiffs(i) Then
            freq = freq + 1
        End If
    Next j
    
    Rfreq(i + 1, 1) = absDiffs(i)
    Rfreq(i + 1, 2) = freq
    
    If freq > maxRankFreq Then
        maxRankFreq = freq
    End If
    
    freq = 0
Next i

'now for the ranks and sum of ranks
nD0 = 0
Rsum = 0
Wmin = 0
Rd0 = 0 'for sum of ranks of differences of 0
ReDim r(1 To nr, 1 To 3)
r(1, 1) = Rfreq(1, 1)
r(1, 2) = Rfreq(1, 2)

If Rfreq(1, 2) = 1 Then
    r(1, 3) = 1
    Else
    r(1, 3) = (1 + 1 + Rfreq(1, 2) - 1) / 2
End If

If scores(0) > hypMed Then
    Rsum = Rsum + r(1, 3)
ElseIf scores(0) = hypMed Then
    nD0 = nD0 + 1
    Rd0 = Rd0 + r(1, 3)
Else
    Wmin = Wmin + r(1, 3)
End If


For i = 2 To nr
    r(i, 1) = Rfreq(i, 1)
    r(i, 2) = Rfreq(i, 2)
    
    If Rfreq(i, 2) = 1 Then
        r(i, 3) = i
    ElseIf Rfreq(i, 1) <> Rfreq(i - 1, 1) Then
        r(i, 3) = (i + i + Rfreq(i, 2) - 1) / 2
    Else
        r(i, 3) = r(i - 1, 3)
    End If
        
    If scores(i - 1) > hypMed Then
        Rsum = Rsum + r(i, 3)
    ElseIf scores(i - 1) = hypMed Then
        nD0 = nD0 + 1
        Rd0 = Rd0 + r(i, 3)
    Else
        Wmin = Wmin + r(i, 3)
    End If
    
Next i

If eqMed = "wilcoxon" Or eqMed = "pratt" Or appr = "exact" Then
    w = Rsum
ElseIf eqMed = "zsplit" Then
    w = Rd0 / 2 + Rsum
End If

If out = "w" Then

    res = w

Else

    If appr = "exact" Then
        If maxRankFreq > 1 Then
            res = "ties occur, cannot compute exact method"
        Else
            statistic = WorksheetFunction.Min(w, Wmin)
            pVal = wcdf(statistic, nr) * 2
            res = pVal
        End If

    Else
    
        If eqMed = "zsplit" Then
            nr = n
        End If
        
        f1 = nr + 1
        s2 = nr * f1 * (2 * nr + 1) / 24
        rAvg = nr * f1 / 4
        
        
        If eqMed = "pratt" Then
            'normal approximation adjustment based on Cureton (1967)
            s2 = s2 - nD0 * (nD0 + 1) * (2 * nD0 + 1) / 24
            rAvg = (nr * f1 - nD0 * (nD0 + 1)) / 4
        End If
        
        
        If ties = True Then
            'ties correction
            t = 0
            For i = 1 To nr - 1
                If Rfreq(i, 1) <> Rfreq(i + 1, 1) Then
                    If eqMed = "wilcoxon" Or eqMed = "pratt" Then
                    'exclude those equal to hypothesized median
                        If Rfreq(i, 1) <> 0 Then
                            t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
                        End If
                    Else
                        t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
                    End If
                    
                End If
            Next i
            If Rfreq(nr - 1, 1) = Rfreq(nr, 1) Then
                t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
            End If
            
            s2 = s2 - t
        End If
        
        
        se = Sqr(s2)
        
        If cc = True Then
            num = Abs(w - rAvg) - 0.5
        Else
            num = Abs(w - rAvg)
        End If
            
        If appr = "imant" Then
            tValue = num / Sqr((s2 * nr - (w - rAvg) ^ 2) / (nr - 1))
            df = nr - 1
            statistic = tValue
            
        Else:
            zValue = num / se
            statistic = zValue
        End If
        
        If appr = "imanz" Then
            zValue = zValue / 2 * (1 + Sqr((nr - 1) / (nr - zValue ^ 2)))
            statistic = zValue
        End If
    End If
    
    If out = "statistic" Then
        res = statistic
    Else
        If appr = "imant" Then
            pVal = WorksheetFunction.T_Dist_2T(Abs(tValue), df)
            res = pVal
        ElseIf appr <> "exact" Then
            pVal = 2 * (1 - WorksheetFunction.Norm_S_Dist(Abs(zValue), True))
            res = pVal
        End If
                
    End If
End If

ts_wilcoxon_os = res

End Function

Function ts_wilcoxon_os_arr(data As Range, _
                    Optional hypMed = "none", _
                    Optional ties = True, _
                    Optional appr = "wilcoxon", _
                    Optional eqMed = "wilcoxon", _
                    Optional cc = False)

testUsed = "one-sample Wilcoxon signed rank test"

If hypMed = "none" Then
    hypMed = (WorksheetFunction.Min(data) + WorksheetFunction.Max(data)) / 2
End If

Dim n, nr As Integer
n = WorksheetFunction.Count(data)

If eqMed = "wilcoxon" Or appr = "exact" Then
    nr = n - WorksheetFunction.CountIf(data, hypMed)
Else
    nr = n
End If

Dim absDiffs() As Double
ReDim absDiffs(0 To nr - 1)
Dim scores() As Double
ReDim scores(0 To nr - 1)
Dim i, k As Integer
k = 0
i = 1
Do While k < nr
    If WorksheetFunction.IsNumber(data(i)) Then
        
        If eqMed = "wilcoxon" Then
            If data(i) <> hypMed Then
                absDiffs(k) = Abs(data(i) - hypMed)
                scores(k) = data(i)
                k = k + 1
            End If
        
        Else
            absDiffs(k) = Abs(data(i) - hypMed)
            scores(k) = data(i)
            k = k + 1
        End If
    End If
    
    i = i + 1
Loop
    

'sort scores based on absolute differences
changes = 1
Do While changes <> 0
    changes = 0
    For i = 1 To nr - 1
        If absDiffs(i - 1) > absDiffs(i) Then
            ff1 = absDiffs(i)
            ff2 = scores(i)
            absDiffs(i) = absDiffs(i - 1)
            scores(i) = scores(i - 1)
            absDiffs(i - 1) = ff1
            scores(i - 1) = ff2
            
            changes = 1
        End If
    Next i
Loop

'we need the ranks for which we need the rank frequencies
'store for each score how often it occurs
'also check if ties actually occur
maxRankFreq = 1
Dim Rfreq As Variant
ReDim Rfreq(1 To nr, 1 To 3)
For i = 0 To nr - 1
    For j = 0 To nr - 1
        If absDiffs(j) = absDiffs(i) Then
            freq = freq + 1
        End If
    Next j
    
    Rfreq(i + 1, 1) = absDiffs(i)
    Rfreq(i + 1, 2) = freq
    
    If freq > maxRankFreq Then
        maxRankFreq = freq
    End If
    
    freq = 0
Next i

'now for the ranks and sum of ranks
nD0 = 0
Rsum = 0
Wmin = 0
Rd0 = 0 'for sum of ranks of differences of 0
ReDim r(1 To nr, 1 To 3)
r(1, 1) = Rfreq(1, 1)
r(1, 2) = Rfreq(1, 2)

If Rfreq(1, 2) = 1 Then
    r(1, 3) = 1
    Else
    r(1, 3) = (1 + 1 + Rfreq(1, 2) - 1) / 2
End If

If scores(0) > hypMed Then
    Rsum = Rsum + r(1, 3)
ElseIf scores(0) = hypMed Then
    nD0 = nD0 + 1
    Rd0 = Rd0 + r(1, 3)
Else
    Wmin = Wmin + r(1, 3)
End If


For i = 2 To nr
    r(i, 1) = Rfreq(i, 1)
    r(i, 2) = Rfreq(i, 2)
    
    If Rfreq(i, 2) = 1 Then
        r(i, 3) = i
    ElseIf Rfreq(i, 1) <> Rfreq(i - 1, 1) Then
        r(i, 3) = (i + i + Rfreq(i, 2) - 1) / 2
    Else
        r(i, 3) = r(i - 1, 3)
    End If
        
    If scores(i - 1) > hypMed Then
        Rsum = Rsum + r(i, 3)
    ElseIf scores(i - 1) = hypMed Then
        nD0 = nD0 + 1
        Rd0 = Rd0 + r(i, 3)
    Else
        Wmin = Wmin + r(i, 3)
    End If
    
Next i

If eqMed = "wilcoxon" Or eqMed = "pratt" Or appr = "exact" Then
    w = Rsum
ElseIf eqMed = "zsplit" Then
    testUsed = testUsed + ", z-split method for equal to hyp. med."
    w = Rd0 / 2 + Rsum
End If

If eqMed = "zsplit" Then
    nr = n
End If

f1 = nr + 1
s2 = nr * f1 * (2 * nr + 1) / 24
rAvg = nr * f1 / 4


If eqMed = "pratt" Then
    testUsed = testUsed + ", Pratt method for equal to hyp. med. (inc. Cureton adjustment for normal approximation)"
    'normal approximation adjustment based on Cureton (1967)
    s2 = s2 - nD0 * (nD0 + 1) * (2 * nD0 + 1) / 24
    rAvg = (nr * f1 - nD0 * (nD0 + 1)) / 4
End If

If ties = True Then
    testUsed = testUsed + ", ties correction applied"
    'ties correction
    t = 0
    For i = 1 To nr - 1
        If Rfreq(i, 1) <> Rfreq(i + 1, 1) Then
            If eqMed = "wilcoxon" Or eqMed = "pratt" Then
            'exclude those equal to hypothesized median
                If Rfreq(i, 1) <> 0 Then
                    t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
                End If
            Else
                t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
            End If
            
        End If
    Next i
    If Rfreq(nr - 1, 1) = Rfreq(nr, 1) Then
        t = t + (Rfreq(i, 2) ^ 3 - Rfreq(i, 2)) / 48
    End If
    
    s2 = s2 - t
End If

If appr = "exact" Then
        If maxRankFreq > 1 Then
            testUsed = "ties occur, cannot compute exact method"
            pVal = "n.a."
            df = "n.a."
            statistic = "n.a."
        Else
            statistic = WorksheetFunction.Min(w, Wmin)
            pVal = wcdf(statistic, nr) * 2
            testUsed = "one-sample Wilcoxon signed rank exact test"
            df = "n.a."
        End If
Else
    se = Sqr(s2)
    
    If cc = True Then
        num = Abs(w - rAvg) - 0.5
    Else
        num = Abs(w - rAvg)
    End If
        
    If appr = "imant" Then
        testUsed = testUsed + ", using Iman's t approximation"
        tValue = num / Sqr((s2 * nr - (w - rAvg) ^ 2) / (nr - 1))
        df = nr - 1
        statistic = tValue
    Else:
        zValue = num / se
        statistic = zValue
        df = "n.a."
    End If
    
    If appr = "imanz" Then
        testUsed = testUsed + ", using Iman's z approximation"
        zValue = zValue / 2 * (1 + Sqr((nr - 1) / (nr - zValue ^ 2)))
        statistic = zValue
    End If
    
    If appr = "imant" Then
        pVal = WorksheetFunction.T_Dist_2T(Abs(tValue), df)
    Else
        pVal = 2 * (1 - WorksheetFunction.Norm_S_Dist(Abs(zValue), True))
    End If
End If

'Results
Dim res(1 To 5, 1 To 5)
res(1, 1) = "W"
res(1, 2) = "statistic"
res(1, 3) = "df"
res(1, 4) = "p-value"
res(1, 5) = "test"
res(2, 1) = w
res(2, 2) = statistic
res(2, 3) = df
res(2, 4) = pVal
res(2, 5) = testUsed

ts_wilcoxon_os_arr = res

End Function


Function srf(k, n)
'signed ranks frequencies

If k < 0 Then
    srf = 0
ElseIf k > WorksheetFunction.Combin(n + 1, 2) Then
    srf = 0

ElseIf n = 1 And (k = 0 Or k = 1) Then
    srf = 1
Else
    srf = srf(k - n, n - 1) + srf(k, n - 1)
End If

End Function

Function wpmf(k, n)
'Wilcoxon Signed Ranks pmf
    wpmf = srf(k, n) / (2 ^ n)
End Function


Function wcdf(k, n)
'Wilcoxon Signed Ranks cdf

Dim i As Integer
Dim p As Double


For i = 0 To k
    p = p + wpmf(i, n)
Next i

wcdf = p

End Function