Estimating Variance as a Function of Treatment Rank Class

Imagine that we have a treatment that we give to five different groups of individuals.  Each individual has a variable response which as a unique mean and variance based on the treatment.  We do not know how the means will change but we believe the variance of responses will expand depending upon what level of treatment the individual gets.  We would like to expressly model both the differences in means and that of the variances.

This code was formulated in response to a question posted on CrossValidated.

We want to solve
$$    \max_{\bf {\hat\beta,\hat\gamma}} (\sum_{i=1}(ln(D(x_i, \hat\mu, \hat\gamma_0+\hat\gamma_1 rank))) $$

# Specify how many individuals are in each of our groups
nobs.group <- 500
 
# Simulate our data
grp1 <- data.frame(values=rnorm(nobs.group,5,1), grp=1)
grp2 <- data.frame(values=rnorm(nobs.group,3,2), grp=2)
grp3 <- data.frame(values=rnorm(nobs.group,6,3), grp=3)
grp4 <- data.frame(values=rnorm(nobs.group,5,4), grp=4)
grp5 <- data.frame(values=rnorm(nobs.group,1,5), grp=5)
 
# Group our data into a single object
mydata <- rbind(grp1,grp2,grp3,grp4,grp5)
 
# Speficy the function to maximize (minimize)
lnp <- function(gamma, x, rank) 
  # I include a negative here because the default option with optim is minimize 
  -sum(log(dnorm(x,gamma[1]*(rank==1)+ 
                   gamma[2]*(rank==2)+
                   gamma[3]*(rank==3)+
                   gamma[4]*(rank==4)+
                   gamma[5]*(rank==5), 
                 gamma[6]+gamma[7]*rank)))
 
ans <- optim(c(
  # Specify initial values for parameters to be estimated
  beta1=1,beta2=1,beta3=1,beta4=1, beta5=1, 
  gamma1=1,gamma2=1), 
  # Specify the function to minimize (maximize)
  lnp, 
  # Input dependent variable as x and the explanatory variable as rank
  x=mydata$values, rank=mydata$grp, 
  # Be sure to inlcude the hessian in the return for 
  # calculating standard errors
  hessian=T)
 
# The standard erros can be estimated using the hessian
stand.error <- sqrt(diag(solve(ans$hessian)))
 
# This will create a nice table of results
cbind(par.est=ans$par, 
      stand.error,
      tstat=ans$par/stand.error,
      pvalue=1-pt(ans$par/stand.error, nrow(mydata)-length(ans$par)),
      CILower=ans$par+stand.error*qt(.05,nrow(mydata)-length(ans$par)),
      CIUpper=ans$par+stand.error*qt(.95,nrow(mydata)-length(ans$par))      
      )
 

 
          par.est stand.error      tstat    pvalue    CILower     CIUpper
beta1   4.9894067  0.04038367 123.550112 0.0000000  4.9229567  5.05585658
beta2   2.0009055  0.10955198  18.264440 0.0000000  1.8206415  2.18116942
beta3   3.7640531  0.19407912  19.394427 0.0000000  3.4447027  4.08340355
beta4   6.4818562  0.23879420  27.144111 0.0000000  6.0889287  6.87478375
beta5  -0.5547626  0.29730735  -1.865957 0.9689176 -1.0439715 -0.06555377
gamma1 -0.4849449  0.05684407  -8.531142 1.0000000 -0.5784798 -0.39140993
gamma2  1.3867000  0.04520519  30.675682 0.0000000  1.3123164  1.46108352