Hi, although the post is 5 years old, this topic is of interest to me. Were you able to get a resolution on the discrepancies between JMP and remlscore? Furthermore JMP gives prediction/confidence intervals whereas remlscore does not. The technical information supplied by JMP is very sparse so any insight would be greatly appreciated.
Here are some useful references:
Smyth, G. K. (2002). “An efficient algorithm for REML in heteroscedastic regression”, Journal of Computational and Graphical Statistics, 11, pp. 1-12.
Verbyla, A.P. (1993). “Modelling Variance Heterogeneity: Residual Maximum Likelihood and Diagnostics”, Journal of the Royal Statistical Society, Series B (Methodological) 55 (2), pp. 493-508.
Hi JMPers!
I'm trying to get a deeper understanding of the loglinear variance models. Therefore I used the InjectionMolding-example from the JMP-documentation. Working with JMP I am able to reproduce the results from the website, of course.
Additionally I implemented an ML-estimator (based on the Harvey-paper) in R. Further I tried the remlscore function from the statmod-package in R (REML estimation like Smyth 2002 describes it). The funny result is, that I get quite similar results with both functions in R for the estimated residual variance (~1.88, ~1.9) but that number is pretty different in JMP: ~6.69.
Now I was curious if this is just some difference in the parametrization? All other numbers match quite nicely.
Here is the R-code, if it helps:
library(statmod)
# Sample Data
d = read.csv("Desktop/InjectionMolding.csv")
x = cbind(rep(1, nrow(d)), d[,c(2, 3)], d[,2]*d[,3]) # MoldTemp, Screw.Speed, 2FI
x = apply(x, 2, as.numeric)
z = d[,4] # Hold.Time
z = cbind(rep(1, length(z)), z)
y = d[,9] # Shrinkage
y = as.matrix(y)
colnames(x) = c("Const", "MoldTemp", "ScrewSpeed", "MoldXScrew")
colnames(z) = c("Const", "HoldTime")
colnames(y) = c("Strength")
reml_est = remlscore(y=y,X=x,Z=z)
reml_est
