R Under development (unstable) (2025-05-10 r88193 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # 
> # test response models
> # 
> 
> require(depmixS4)
Loading required package: depmixS4
Loading required package: nnet
Loading required package: MASS
Loading required package: Rsolnp
Loading required package: nlme
> 
> # binomial response model
> x <- rnorm(1000)
> # library(boot)
> invlogit <- function(lp) {exp(lp)/(1+exp(lp))}
> p <- invlogit(x)
> ss <- rbinom(1000,1,p)
> mod <- GLMresponse(cbind(ss,1-ss)~x,family=binomial())
> fit(mod)
Model of type binomial (logit), formula: cbind(ss, 1 - ss) ~ x
Coefficients: 
(Intercept)           x 
  0.1572014   1.0534315 
Probality at zero values of the covariates. 
0.5392196 
> glm(cbind(ss,1-ss)~x, family=binomial)

Call:  glm(formula = cbind(ss, 1 - ss) ~ x, family = binomial)

Coefficients:
(Intercept)            x  
     0.1572       1.0534  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:	    1385 
Residual Deviance: 1180 	AIC: 1184
> 
> # gamma response model
> x=runif(1000,1,5)
> res <- rgamma(1000,x)
> ## note that gamma needs proper starting values which are not
> ## provided by depmixS4 (even with them, this may produce warnings)
> mod <- GLMresponse(res~x,family=Gamma(),pst=c(0.8,1/0.8))
> fit(mod)
Model of type Gamma (inverse), formula: res ~ x
Coefficients: 
(Intercept)           x 
  0.7486414  -0.1179044 
There were 17 warnings (use warnings() to see them)
> glm(res~x,family=Gamma)

Call:  glm(formula = res ~ x, family = Gamma)

Coefficients:
(Intercept)            x  
     0.7486      -0.1179  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:	    579.6 
Residual Deviance: 425 	AIC: 3641
> 
> # multinomial response model
> x <- sample(0:1,1000,rep=TRUE)
> mod <- GLMresponse(sample(1:3,1000,rep=TRUE)~x,family=multinomial(),pstart=c(0.33,0.33,0.33,0,0,1))
> mod@y <- simulate(mod)
> fit(mod)
Model of type multinomial (mlogit), formula: sample(1:3, 1000, rep = TRUE) ~ x
Coefficients: 
            1           2           3
(Intercept) 0 -0.06771816 -0.01801858
x           0 -0.02140969  0.73035989
Probalities at zero values of the covariates.
0.3428571 0.3204082 0.3367347 
> colSums(mod@y[which(x==0),])/length(which(x==0))
[1] 0.3428571 0.3204082 0.3367347
> colSums(mod@y[which(x==1),])/length(which(x==1))
[1] 0.2529412 0.2313725 0.5156863
> 
> # multivariate normal response model
> mn <- c(1,2,3)
> sig <- matrix(c(1,.5,0,.5,1,0,0,0,2),3,3)
> y <- mvrnorm(1000,mn,sig)
> mod <- MVNresponse(y~1)
> fit(mod)
Multivariate Normal Model, formula: y ~ 1
Coefficients: 
                 [,1]     [,2]     [,3]
(Intercept) 0.9821976 1.975298 3.027757
Sigma: 
            [,1]       [,2]        [,3]
[1,]  0.97133045 0.52179815 -0.04602515
[2,]  0.52179815 1.04265511  0.00819608
[3,] -0.04602515 0.00819608  1.95599781
> colMeans(y)
[1] 0.9821976 1.9752979 3.0277572
> var(y)
            [,1]        [,2]         [,3]
[1,]  0.97230275 0.522320467 -0.046071219
[2,]  0.52232047 1.043698804  0.008204284
[3,] -0.04607122 0.008204284  1.957955765
> 
> # normal (gaussian) response model
> y <- rnorm(1000)
> mod <- GLMresponse(y~1)
> fm <- fit(mod)
> cat("Test gaussian fit: ", all.equal(getpars(fm),c(mean(y),sd(y)),check=FALSE))
Test gaussian fit:  names for target but not for current Mean relative difference: 0.0004916515> 
> 
> # poisson response model
> x <- abs(rnorm(1000,2))
> res <- rpois(1000,x)
> mod <- GLMresponse(res~x,family=poisson())
> fit(mod)
Model of type poisson (log), formula: res ~ x
Coefficients: 
(Intercept)           x 
 -0.3926746   0.4784137 
> glm(res~x, family=poisson)

Call:  glm(formula = res ~ x, family = poisson)

Coefficients:
(Intercept)            x  
    -0.3927       0.4784  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:	    1585 
Residual Deviance: 1127 	AIC: 3240
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.87    0.15    1.01