I am creating an R code from scratch that is replicating an existing model in VBA that loops over improvements using Newton Raphson several times until the output differences stabilise. As such any small rounding differences on what is being inputted into the start of the code end up causing it to fall over.

This is currently happening in that the inputs E are being brought in through read.excel on the original dataset in excel that the VBA uses. The values are for example 365318.425165619 and 313266.786807882 where they are displayed as 9 decimal places. Then in the R dataset that has been read in from these they are giving the values as 365318.42516561900266 and 313266.78680788201746 where they are displayed as 14 decimal places.

My question is that whether there is something in R that operates differently to excel and whether this is causing my calculation to fall over further down the calculations as the initial calculations seem to be fine.

For reference it is then used here in the alpha calculation that then gives slightly different output in the R dataset (-3.98568702613257741518, -3.98577596443294002171) compared to the VBA (-3.98568702613256, -3.98577596443293)

Alpha R Calculation Steps:

#creating dobj for alpha deviance
APCI_summary$dobj <- 2 * APCI_summary$weight * (APCI_summary$final_exposure * APCI_summary$m_xy - APCI_summary$calc_deaths)
dalphaobj_values <- aggregate(APCI_summary$dobj, by = list(APCI_summary$age), FUN = sum)$x
dalphaobj_m <- matrix(dalphaobj_values, nrow = length(dalphaobj_values), ncol = 1)

#creating d2obj for alpha deviance
APCI_summary$d2obj <- 2 * APCI_summary$weight * APCI_summary$final_exposure * APCI_summary$m_xy
d2alphaobj_values <- aggregate(APCI_summary$d2obj, by = list(APCI_summary$age), FUN = sum)$x
d2alphaobj_m <- diag(d2alphaobj_values)

#creating dobj and d2obj for alpha penalty
alpha_table <- data.frame(age = seq(age_min, age_max), alpha_x = APCI_summary$alpha_x[APCI_summary$year == year_max])
alpha_x_values <- alpha_table$alpha_x

# Convert the extracted column into a matrix
Alpha <- matrix(alpha_x_values, nrow = length(alpha_x_values), ncol = 1)
tAlpha <- t(Alpha)

dP_temp_alpha <- 2 * lambda_alpha * tDiff_3 %*% Diff_3 %*% Alpha

d2P_temp_alpha <- 2 * lambda_alpha * tDiff_3 %*% Diff_3

#combining dobj and d2obj respectively
dalpha_total_m <- dalphaobj_m + dP_temp_alpha
d2alpha_total_m <- d2alphaobj_m + d2P_temp_alpha

#final update step including the update to alpha values
d2Obj_inv <- solve(d2alpha_total_m) %*% dalpha_total_m
Alpha <- Alpha - d2Obj_inv
Alpha_data <- data.frame(age = seq(age_min, age_max), alpha_x = Alpha)

#then needs to update the calculation for this new Alpha value
APCI_summary <- APCI_summary %>%
  left_join(Alpha_data, by = "age") %>%
  mutate(alpha_x = coalesce(alpha_x.y, alpha_x.x)) %>%
  select(-alpha_x.x, -alpha_x.y)

APCI_summary$log_m_xy <- log_m_xy(APCI_summary$alpha_x, APCI_summary$beta_x, APCI_summary$year, APCI_summary$kappa_x, APCI_summary$gamma_x)
APCI_summary$m_xy <- exp(APCI_summary$log_m_xy)

APCI_summary$deviance <- APCI_summary$weight * (APCI_summary$calc_deaths * log(APCI_summary$calc_deaths)
                                                - APCI_summary$calc_deaths
                                                - APCI_summary$calc_deaths * log(APCI_summary$final_exposure * APCI_summary$m_xy)
                                                + APCI_summary$final_exposure * APCI_summary$m_xy)

Edit: this is also the VBA code it is trying to replicate which could possibly flag why it is falling over itself:

Private Sub ImproveAlpha(Lambda)

Dim X, Y, W, C, dlogm, X1, X2
Dim dObj(), d2Obj() As Double    ' first and second differentials of the Objective function
ReDim dObj(RP.MinAge To RP.MaxAge, 1 To 1)
ReDim d2Obj(RP.MinAge To RP.MaxAge, RP.MinAge To RP.MaxAge) As Double

' contribution from deviance
For X = RP.MinAge To RP.MaxAge
    For Y = RP.MinYear To RP.MaxYear
        W = RP.Weight(Y)
        C = Y - X
        If RP.MinCohort <= C And C <= RP.MaxCohort Then
            dlogm = 1
            dObj(X, 1) = dObj(X, 1) + 2 * W * (E(X, Y) * m(X, Y) - D(X, Y)) * dlogm
            d2Obj(X, X) = d2Obj(X, X) + 2 * W * (E(X, Y) * m(X, Y)) * dlogm ^ 2
        End If
    Next Y
Next X


' contribution from regularisation penalties
Dim Diff, P, dP, d2P, dP_temp

Diff = DifferenceMatrix(RP.NumAges, ORDERALPHA)
P = MMult(Transpose(Diff), Diff)

ReDim dP(RP.MinAge To RP.MaxAge, 1 To 1)
ReDim d2P(RP.MinAge To RP.MaxAge, RP.MinAge To RP.MaxAge)
dP_temp = MMult(P, VectorToMatrix(Alpha))


For X1 = RP.MinAge To RP.MaxAge
    dObj(X1, 1) = dObj(X1, 1) + 2 * Lambda * dP_temp(X1 + 1 - RP.MinAge, 1)
    For X2 = RP.MinAge To RP.MaxAge
        d2Obj(X1, X2) = d2Obj(X1, X2) + 2 * Lambda * P(X1 + 1 - RP.MinAge, X2 + 1 - RP.MinAge)
    Next X2
Next X1


' update the parameters
Dim Delta
Delta = MMult(MInverse(d2Obj), dObj)
For X = RP.MinAge To RP.MaxAge
    Alpha(X) = Alpha(X) - Delta(1 + X - RP.MinAge, 1)
Next X

End Sub