[R] Help me improving my code

[David Winsemius]

This is rather obviously homework, and you have not read the Posting  
Guide, and you have not addressed the question of academic integrity  
policies that are probably in force at your university.

-- 
David
On Oct 31, 2009, at 1:30 PM, md. jakir hussain talukdar wrote:

> Hi,
>
> I am new to R. My problem is with the ordered logistic model. Here  
> is my
> question:
>
>
> Generate an order discrete variable using the variable
>
> wrwage1 = wages in first full calendar quarter after benefit  
> application
>
> in the following way:
> *
>
> wage*1*Ordered *=
>
> 1 *if*0 *· wrwage*1 *< *1000
>
> 2 *if*1000 *· wrwage*1 *< *2000
>
> 3 *if*2000 *· wrwage*1 *< *3000
>
> 4 *if*3000 *· wrwage*1 *< *4000
>
> 5 *ifwrwage*1 *¸ *4000.
>
> Estimate an order logit model using the new created discrete  
> variable as the
>
> dependent variable with the following control variables:
>
> black = 1 if black
>
> hispanic = 1 if hispanic
>
> otherrace = 1 if black = 0 and hispanic = 0
>
> female = 1 if female
>
> wba = claimant's UI weekly bene¯t amount (*=week*)
>
> age = age of claimant when claim is ¯led (years)
>
> bp1000 = base period wages (thousand $)
>
> inuidur1 = duration of unemployment (weeks)
>
> I have written the program as:
>
>> pennw <- read.table("G:/Econometrics II/pennsylvania.txt",  
>> header=TRUE)
>
>> names(pennw)
>
> [1] "group"     "t0"        "t1"        "t2"        "t3"        "t4"
>
> [7] "intercept" "black"     "hispanic"  "female"    "wba"        
> "bonus"
>
> [13] "age"       "recall"    "pd"        "bp1000"    "numpay"     
> "wrwage0"
>
> [19] "wrwage1"   "wrwage2"   "wrwage3"   "wrwage4"   "inuidur1"   
> "regui"
>
>> wage <- pennw$wrwage1
>
>> wage[wage>=0 & wage<1000] <- 1
>
>> wage[wage>=1000 & wage<2000] <- 2
>
>> wage[wage>=2000 & wage<3000] <- 3
>
>> wage[wage>=3000 & wage<4000] <- 4
>
>> wage[wage>=4000] <- 5
>
>> intercept <- pennw$intercept
>
>> blac <- pennw$black
>
>> hisp <- pennw$hispanic
>
>> fem <- pennw$female
>
>> wba <- pennw$wba
>
>> age <- pennw$age
>
>> bpw <- pennw$bp1000
>
>> dur <- pennw$inuidur1
>
>> fw <- factor(wage)
>
>> Ik <- model.matrix(~fw-1)
>
>> I1 <- Ik[,1]
>
>> I2 <- Ik[,2]
>
>> I3 <- Ik[,3]
>
>> I4 <- Ik[,4]
>
>> I5 <- Ik[,5]
>
>> yelmon <- function(theta,x){
>
> + mu1 <- theta[1]
>
> + mu2 <- theta[2]
>
> + mu3 <- theta[3]
>
> + mu4 <- theta[4]
>
> + mu5 <- theta[5]
>
> + beta <- theta[6]
>
> + logL <-
> sum((Ik*log((exp(mu2-(x*beta))/(1+exp(mu2-(x*beta))))-(exp(mu1- 
> (x*beta))/(1+exp(mu1-(x*beta))))))+(Ik*log((exp(mu3-(x*beta))/ 
> (1+exp(mu3-(x*beta))))-(exp(mu2-(x*beta))/(1+exp(mu2-(x*beta))))))+ 
> (Ik*log((exp(mu4-(x*beta))/(1+exp(mu4-(x*beta))))-(exp(mu3-(x*beta))/ 
> (1+exp(mu3-(x*beta))))))+(Ik*log((exp(mu5-(x*beta))/(1+exp(mu5- 
> (x*beta))))-(exp(mu4-(x*beta))/(1+exp(mu4-(x*beta))))))+(Ik*log(1- 
> (exp(mu5-(x*beta))/(1+exp(mu5-(x*beta)))))))
>
> + return(-logL)
>
> + }
>
>> optim(yelmon, c(0,1), x=c(intercept, blac, hisp, fem, wba, age,  
>> bpw, dur))
>
> I should be getting the values for all beta’s and mu’s here
> ) The program
> should be done by maximum likelihood not by calling ready-made  
> commands
>
> Rubel
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
Heritage Laboratories
West Hartford, CT