Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I have following data stored in a file. I am applying 'glm' in R to find linear regression equation to best predict the 'output'.

> tmpData
   logOfOutput randomSample multiplied part1 part2 randNormalMean100Std20 output
1    0.0000000           33         11     1    19               89.65387      1
2    0.6931472           76         24     2    18              128.23471      2
3    1.0986123           12         39     3    17              103.70930      3
4    1.3862944           68         56     4    16               99.12617      4
5    1.6094379           50         75     5    15               95.68173      5
6    1.7917595            7         96     6    14              129.27551      6
7    1.9459101           70        119     7    13              104.59333      7
8    2.0794415           55        144     8    12              102.15247      8
9    2.1972246           20        171     9    11               72.43795      9
10   2.3025851           24        200    10    10               80.63634     10
11   2.3978953           32        231     9    11              105.03423     11
12   2.4849067           97        264     8    12               78.10613     12
13   2.5649494           28        299     7    13              107.95286     13
14   2.6390573           99        336     6    14               80.07396     14
15   2.7080502           66        375     5    15              102.01156     15
16   2.7725887           95        416     4    16              119.07361     16
17   2.8332133           42        459     3    17               64.19354     17
18   2.8903718           53        504     2    18              106.23402     18
19   2.9444390           85        551     1    19              151.07976     19
20   2.9957323           48        600     0    20               82.78324     20

I am using the following code to perform the same

fn = "delnowSample.txt"
tmpData = read.table(fn, header = TRUE,  sep= "\t" , blank.lines.skip = TRUE)
cnames = colnames(tmpData)
(fmla <- as.formula(paste(cnames[length(cnames)], " ~ ", paste(cnames[1:(length(cnames)-1)],collapse= "+")))  )
model <- try(glm(formula = fmla, family=binomial(), na.action=na.omit, data=tmpData));
summary(model) 

The output that I get is as follow:

> summary(model)

Call:
glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse = "+"))), 
    family = gaussian(), na.action = na.omit)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.37926  -0.11242  -0.03441   0.16087   0.28200  

Coefficients: (1 not defined because of singularities)
                                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                                0.2638036  0.3078536   0.857  0.40592    
unlist(tmpData["logOfOutput"])             0.9202273  0.2727884   3.373  0.00455 ** 
unlist(tmpData["randomSample"])            0.0026201  0.0018177   1.441  0.17145    
unlist(tmpData["multiplied"])              0.0288073  0.0012359  23.308 1.34e-12 ***
unlist(tmpData["part1"])                   0.2106002  0.0403442   5.220  0.00013 ***
unlist(tmpData["part2"])                          NA         NA      NA       NA    
unlist(tmpData["randNormalMean100Std20"]) -0.0006214  0.0024922  -0.249  0.80673    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for gaussian family taken to be 0.04403284)

    Null deviance: 665.00000  on 19  degrees of freedom
Residual deviance:   0.61646  on 14  degrees of freedom
AIC: 1.1676

Number of Fisher Scoring iterations: 2

To a large extent it is predicting the Pr(z) correctly as we can see the probabilities of random variable are not significant. The R-square is also high (1-residual.deviance/null.deviance), close to 1.

Question 1: In the above data 'part1+part2' is equal to output variable. Is 'glm' not able to identify such type of relations?

Question 2: Why the degree of freedom of null and residual deviance are different?

Question 3: I need to convert the output variable into categorical variable (i.e. Everything <=10 is 'no' and more than this is 'yes'). What is the best way to call 'glm', when the response variable is 'categorical'. I tried converting 'no' to '0' and 'yes' to 1, and called glm as follows:

 model <- try(glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse= "+"))), family=binomial(), na.action=na.omit));

I am getting warning message with this code. Also, I am not sure if this is the correct way to call categorical variable.

Edit:

I have the following categorical data:

> tmpData
           x1  x2 x3  y1
1  0.16294456   1  1  no
2  0.80494934   2  2  no
3  0.28962222   1  3  no
4  0.07177347   2  4  no
5  0.54830544   1  5  no
6  0.67655327   2  6  no
7  0.45189608   1  7  no
8  0.82412502   2  8  no
9  0.09076793   1  9  no
10 0.12221227   2 10  no
11 0.56751754 111 11 yes
12 0.04970992 222 12 yes
13 0.56162037 111 13 yes
14 0.96617891 222 14 yes
15 0.50994534 112 15 yes
16 0.70093692 212 16 yes
17 0.02034940 212 17 yes
18 0.78356903 121 18 yes
19 0.58439662 213 19 yes
20 0.31729282 212 20 yes

And the following code:

  fn = "delnowSample.txt"
  tmpData = read.table(fn, header = TRUE,  sep= "\t" , blank.lines.skip = TRUE)
  tmpData
  model <- glm(formula = 'y1~x1+x2+x3', family=binomial(), na.action=na.omit, data=tmpData)
  summary(model) 

This one doesn't seem to be working??

share|improve this question
    
What's the unlisting for? It makes it very unclear to see what the code is doing (and therefore to answer your question). I suggest constructing the variables you need in a new data frame and calling them directly. –  Aaron Jan 8 '13 at 20:55
    
@Aaron, thanks. I am using 'unlist' as glm was considering these columns as list and was showing me "invalid type (list) for variable 'tmpData["output"]" error and that is why I used unlist. –  user1140126 Jan 8 '13 at 21:17
    
You should not need the unlist after read.table() because it should return a dataframe. If you do str() on your tmpData, it should give you more info on what happen at that stage. You could also bypass completely the call to as.formula() by calling glm([your formula], data=tmpData). It is still possible to build [your formula] with paste() and the colnames if you don't want to type it in. That should save you some of the trouble and make the code easier to read. –  Antoine Vernet Jan 8 '13 at 21:48
1  
Yes to what @AntoineVernet says, or if you do need to make the formula automatically, xn <- cnames[1:(length(cnames)-1)] (and similarly for dep) should be fine, if you add the data=tmpData parameter. (Or just use ``[[` instead of [...) –  Aaron Jan 8 '13 at 21:55
    
Also, use lm instead of glm with the normal error. –  Aaron Jan 8 '13 at 21:56
show 3 more comments

1 Answer

up vote 4 down vote accepted
  • In response to question 1

part1 + part2 =20 not output. R has recognized this.

It states

Coefficients: (1 not defined because of singularities)

and does not attempt to define a coefficient for part2 (Hence the NA values)

unlist(tmpData["part2"])                          NA         NA      NA       NA    

That being said, your code is not consistent as you define

model <- try(glm(formula = fmla, family=binomial(), na.action=na.omit, data=tmpData));

Which should give an error, if fmla is

output ~ logOfOutput +randomSample+ multiplied part1 +part2+ randNormalMean100Std20

as your code suggests as output appears to be a continuous variable, to quote from ?glm

For binomial and quasibinomial families the response can also be specified as a factor (when the first level denotes failure and all others success) or as a two-column matrix with the columns giving the numbers of successes and failures.

In addition why are you using logOfOutput in a linear model to predict output?

Then you present the output of a completely different model (gaussian error, different formula as given away by the Call and the awful mess of coefficient names)

summary(model)


Call:
glm(formula = as.formula(paste(dep, " ~ ", paste(xn, collapse = "+"))), 
    family = gaussian(), na.action = na.omit
  • In response to question 2.

The degrees of freedom are different as the null deviance is the deviance for a model contains only the intercept term.

The residual deviance is the deviance from the model you actually fitted (with 4 more parameters being estimated).

  • Question 3

Although I don't think there is much point going much further until you can create a small reproducible example and you are clear on what your data are and your model should be.

The issue with your second example data set appears data related.

share|improve this answer
    
thanks for your elaborate and well explained answer. I admit that I messed up the code while editing the post. –  user1140126 Jan 9 '13 at 15:02
    
I posted the same question "stats.stackexchange.com/questions/47327/…; –  user1140126 Jan 9 '13 at 15:38
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.