I am not very experienced in converting equations to code. I got stuck on converting the score function of the partial likelihood to code in julia to be evaluated in JuMP.
The score function which, when solved for beta at 0, is the max.
I made a small simple data set.
Using DataFrames, DataFramesMeta, JuMP, Ipopt
#build DataFrame
times = [6,7,10,15,19,25]
is_censored = [1,0,1,1,0,1]
x= is_control = [1,1,0,1,0,0]
m = Model(solver=IpoptSolver(print_level=0))
using DataFrames
df = DataFrame();
df[:times]=times;
df[:is_censored]= is_censored;
df[:x]=x;
df
#sort df
df_sorted = sort!(df, cols = [order(:times)])
#make df_risk and df_uncensored
df_uncensored = @where(df_sorted, :is_censored .== 0)
df_risk = df_sorted
#use JuMP
##convert df to array
uncensored = convert(Array,df_uncensored[:x])
risk_set = convert(Array,df_risk[:x])
risk_index = convert(Array,find(is_censored .== 0))
x = convert(Array, x)
@variable(m, β, start = 0.0)
# score
@NLobjective(m, Max, sum(uncensored[i] - (([sum(exp(risk_set[j]*β)*x[j]) for j = risk_index[i]:length(risk_set)]) / ([(sum(exp(risk_set[j]*β)*x[j])) for j=risk_index[i]:length(risk_set)])) for i = 1:length(uncensored)))
The error I am getting is
ERROR: exp is not defined for type AffExpr. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.
Stacktrace:
[1] exp(::JuMP.GenericAffExpr{Float64,JuMP.Variable}) at /home/icarus/.julia/v0.6/JuMP/src/operators.jl:630
[2] collect(::Base.Generator{UnitRange{Int64},##58#60}) at ./array.jl:418
[3] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/parseExpr_staged.jl:489 [inlined]
[4] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/parsenlp.jl:226 [inlined]
[5] macro expansion at /home/icarus/.julia/v0.6/JuMP/src/macros.jl:1086 [inlined]
[6] anonymous at ./<missing>:?
The error says the the exponential is the problem, but I previously did the log likelihood which has an exponential in it and got no errors.
In R the beta = -1.3261
If the above code was working I would expect that same output after running
solve(m)
println("β = ", getvalue(β))
The code in the question is a bit convoluted, but the following is an attempt to extract the relevant parts for a working method:
using JuMP, Ipopt
times = [6,7,10,15,19,25];
is_censored = 1-[1,0,1,1,0,1];
is_control = 1-[1,1,0,1,0,0];
uncensored = find(is_censored .== 0)
println("times = $times")
println("is_censored = $is_censored")
println("is_control = $is_control")
m = Model(solver=IpoptSolver(print_level=0))
@variable(m, β, start = 0.0)
@NLobjective(m, Max, sum(log(1+(-1)^is_control[uncensored[i]]*
sum((-1)^is_control[j]*exp(is_control[j]*β) for j=uncensored[i]:length(times))/
sum( exp(is_control[j]*β) for j=uncensored[i]:length(times)))
for i=1:length(uncensored)))
solve(m)
println("β = ", getvalue(β))
This outputs:
times = [6,7,10,15,19,25]
is_censored = [0,1,0,0,1,0]
is_control = [0,0,1,0,1,1]
β = -1.3261290591982942
The β is the same one from the question, so I guess the tweaks to the input are correct and the formula is the log-likelihood. The beginning of the expression inside the log uses a common trick of choosing a +1 or -1 sign according to a 0/1 value bool with (-1)^bool.