Ciao, anche se il post ha 5 anni, questo argomento mi interessa. Sei riuscito a risolvere le discrepanze tra JMP e remlscore? Inoltre JMP fornisce intervalli di previsione / confidenza mentre remlscore no. Le informazioni tecniche fornite da JMP sono molto scarse, quindi qualsiasi intuizione sarebbe molto apprezzata.
Ecco alcuni riferimenti utili:
Smyth, GK (2002). "Un algoritmo efficiente per REML nella regressione eteroschedastica", Journal of Computational and Graphical Statistics, 11, pp. 1-12.
Verbyla, AP (1993). "Modeling Variance Heterogeneity: Residual Maximum Likelihood and Diagnostics", Journal of the Royal Statistical Society, Series B (Methodological) 55 (2), pp. 493-508.
Ciao JMPers!
Sto cercando di ottenere una comprensione più approfondita dei modelli di varianza loglineare. Pertanto ho utilizzato l'esempio InjectionMolding dalla documentazione di JMP. Naturalmente, lavorando con JMP sono in grado di riprodurre i risultati dal sito web.
Inoltre ho implementato uno stimatore ML (basato sull'Harvey-paper) in R. Inoltre ho provato la funzione remlscore dal pacchetto statmod in R (la stima REML come la descrive Smyth 2002). Il risultato divertente è che ottengo risultati abbastanza simili con entrambe le funzioni in R per la varianza residua stimata (~ 1,88, ~ 1,9) ma quel numero è abbastanza diverso in JMP: ~ 6,69.
Ora ero curioso se questa è solo qualche differenza nella parametrizzazione? Tutti gli altri numeri corrispondono abbastanza bene.
Ecco il codice R, se aiuta:
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
