title: Please help me refine this zero-inflated negative binomial model

text: I have been working on a baseball model to predict success at the major league level using minor league statistics. After posting multiple threads on this site (1, 2, 3) and receiving valuable feedback, I have settled on a zero- inflated negative binomial model as being the best fit for my data. For those who do not want to go back through old threads, I will recap some of the story here. Also, for those who have read the old threads, some of the details regarding the variables I am using have changed. In my model, the dependent variable, offensive career wins above replacement (oWAR), is a proxy for success at the MLB level and is measured as the sum of offensive contributions for every play the player is involved in over the course of his career (details here - http://www.fangraphs.com/library/misc/war/). The independent variables are z-scored minor league offensive variables for statistics that are thought to be important predictors of success at the major league level including age (players with more success at a younger age tend to be better prospects), strike out rate [SOPct], walk rate [BBPct] and adjusted production (a global measure of offensive production). Additionally, since position is an important determinant of whether a players makes the major leagues (those who play at easier positions will be required to perform at a higher offensive level in order to have the same value as a player at a more difficult position), I have included dummy variables to account for position. Note that I have not included the position dummy in the count portion of the model as the oWAR dependent variable has already been adjusted for the difficulty level of the position played by the player. **EDIT: I have added the paragraphs below in response to the following comments in the answer below:** _"I see you do not include the same covariates in the Logit and the negative binomial process - why not? Ususally, each relevant predictor would be expected to influence both processes."_ I think it would help if I explained the data generating process. A player plays in the minor leagues. At some point, when they have demonstrated enough skill at the minor league level (this is a combination of statistical success and observed traits that scouts believe will allow them to be successful at the major league level), a player is promoted to the major leagues. At this point, they have an opportunity to accrue oWAR. At this point, the first data- generating process (captured by the logit model) ends. Now, a different data generating process takes over, whereby players accumulate oWAR depending on how they play at the major league level. Some players will not perform well and accumulate zero oWAR, the same as a player who did not make the majors. That is one of the reasons I think this model is appropriate. It is not necessarily easy to separate a player who accumulates zero because they aren't good enough to make the major leagues from a player who makes the major leagues but does not succeed at that level (and still ends his career with zero oWAR). I have not included the positional dummy in the count part of the model because the oWAR measure is already adjusted for the position played by the player whereas the minor league statistics are not. When I tried testing them in the model, they were, not surprisingly, not significant. I omitted the BB Pct statistic from the logit part of the model as it was not significant (p = 0.22) _"Since your data appears to be a panel (you observe players/teams repeatedly), you can think about more sophisticated stuff like fixed or random effects. "_ In terms of how the data is sampled, a player can be in the dataset once (if they spend only a year at the level of the minor leagues) or multiple times if they take multiple years to advance. After reading up on fixed vs random effects model, I don't see how I can used a fixed model to predict out of sample players. However, I am sure that there are fixed effects (effects determined by the player that are not captured in the dependent variables) so I don't fully understand how to handle that issue. **END EDIT** After trying a linear, Box-Cox transformed and basic GLM model with generally poor results, I was directed to the zero-inflated negative binomial distribution set of models. After trying out different combination of variables and following the steps in this excellent step-by-step guide for regression models for count data in R, I settled on the following model (shown below). **Model** ![model1](http://i.imgur.com/neOP2RE.jpg) Furthermore, when I re-estimate the standard errors using sandwich standard errors, the model still appears to have appropriate independent variables. **Sandwich Standard Errors** ![sandwich](http://i.imgur.com/47JjOae.jpg) At this point, I do not think I will find a better model type given the dataset I have. However, I am still left with some issues. The first issue may be a function of the dataset I am using. In most examples I have seen that discuss zero-inflated data, there are clearly more zeros than other values (hence, the name). However, the number of zeros still appears to be less than 50% of the total dependent variables and usually not even that high. In my dataset, approximately 87% of the dependent variables are zero i.e., it is hard to have success in major league baseball. I am guessing the model should technically be able to account for this scenario (albeit with less predictive value than a model with more non-zeros) but I am not sure how to check if that is the case. When I create a plot of fitted values and Pearson residuals, and a plot of fitted values and raw residuals, they appear as below: **Fitted vs Pearson residuals** ![FvP](http://i.imgur.com/tlx6ibn.jpg) **Fitted vs raw residuals** ![FvR](http://i.imgur.com/vL8OIcV.jpg) Not knowing exactly what these plots look like in a good-fitting regression, I decided to take the sample data described here and examine the plots in an example where I know the fit has been deemed to be good. **Fitted vs Pearson residuals - Sample problem** ![FvPEx](http://i.imgur.com/TJh2zHF.jpg) **Fitted vs raw residuals - Sample problem** ![FvREx](http://i.imgur.com/YJMxJl9.jpg) Clearly, these plots do not look that similar to mine. I am not sure how much has to do with a) model misspecification b) the fact that my dependent variable has 87% zeros and c) the fact that this is a simple sample problem designed to perfectly fit this model whereas my data is messy, real world data. Any thoughts on this issue would be appreciated. My second issue, which I am not sure if I should be tackling after or simultaneously with the first issue, has to do with the functional form specification. I don't know if my independent variables are in the right form. It has been suggested to me by a friend that I could try a) multiple fractional polynomials with loops or b) informally play around with adding polynomials of covariates, interactions, etc. My issue at this point is that I do not know how to implement point a in R and and I am not sure which forms to try for point b besides randomly choosing some. Once again, help on this separate (but related?) issue would be greatly appreciated. If anyone has any questions, I will do my best to answer them. In my first post (1), I mentioned I could not provide the dataset but I have been given permission to do so if anyone wants to take a look. Thanks again.


title: Confused about standardization

text: I read a lot on the web, but I am still not sure whether I understood completely when we standardize the data (so that it is zero mean unit variance). So, let's say that I have a set of genes and expression levels of those genes among different cancer patients, and I want to cluster the genes into pathways using expression data. What I need to do as a preprocessing step is to standardize each gene so that it is zero mean unit variance across samples, right? That's what people do in the papers I read. Okay, now let's say I want to work on the transposed matrix, that's, let's say I want to cluster patients into different cancer subtypes. Then, should I standardize the data this time for each sample so that each sample is zero mean unit variance across genes? Or, should I again standardize for each gene as before when we were clustering the genes? The latter approach seems to make more sense since gene activity levels are actually dependent to each other, but samples are usually independent. Making genes independent from each other using the former step does not seem logical to me. Could you please clarify on what dimension of the data and when I need standardizing? Thanks.


title: How to create a glmm for continuous and categorical explanatory variables and binary response?

text: After weeks of reading and trying I decided to post my question here because I could not find a convincing solution. I radio tracked two animals for several months and now I want to find out 1) what influences the activity of the animals and 2) when (hour after sunset) they are showing the highest activity and 3) what influences the travel distance. For question 1) I created a generalized linear mixed model, with animal as a random factor looking like this: glmer(cbind(active,inactive)~offspring+season+observation_time+temperature+precipitation+season:temperature+season:precipitation+temperature:precipitation,family=binomial) offspring and season are factors coded with 0 and 1,observation_time in minutes, temperature is in °C and precipitation in mm. The first lines of the summary() are showing that: AIC BIC logLik deviance 438.3251 460.0690 -209.1625 418.3251 Random effects: Groups Name Variance Std.Dev. Name (Intercept) 0.06594 0.2568   Number of obs: 65, groups: Name, 2 So my question is: is this model build up correctly? Is there an improvement possible or neccessary? It is very difficult to work with this data, because most diagnostic plots which are usually used to evaluate models are different because of the random factor. I also wanted to boxcox transform the response but one animal showed no activity one day (the activity is zero) and therefore this is not possible. I try to eliminate variables or interaction terms but in most cases I can only eliminate one interaction term. After that all variables and interaction terms seems to be significant. For variable selection I use differend aproaches (AIC, BIC, AVOVA). For the 3) question I created a LMM like this: lmer(sqrt(distance)~aktivity+season+offspring+observation_time+temperature+precipitation+season:temperature+...(other interaction terms) The result from summary(): REML criterion at convergence: 513.0257 Random effects: Groups Name Variance Std.Dev. Name (Intercept) 11.6 3.405   Residual 295.8 17.199   Number of obs: 62, groups: Name, 2 Does any of this values tell me something about the goodness of my model? How can I check if data transformation is neccessary? And if yes, which one? I'm honest, I'm quite new in this field and all I learnd about R and statistics do not really work with glmm or lmm. Or it's to complicated for me. I also created a gam without the random factor to check the relationship between the variables and response but I don't know what to do with the results (seems to be no linear relationship between activity and observation_time and rainfall). How do I fit variables to my model? An other idea is to fit the model without the random factor and add the random factor afterwards. Would it be ok to do it like this? For the second question - at which hour after sunset they are showing the highest activity - I have no idea how the model could be build up... Sorry for the amount of questions but I'm working for weeks on this and it is very frustrating... Thanks in advance for all your ideas and help! Iris


title: Problem calculating, interpreting regsubsets and general questions about model selection procedure

text: I want to select models using `regsubsets()`. I have a dataframe called olympiadaten (data uploaded: http://www.sendspace.com/file/8e27d0). I first attach this dataframe and then start to analyze, my code is:               attach(olympiadaten)          library(leaps)     a<-regsubsets(Gesamt ~ CommunistSocialist + CountrySize + GNI + Lifeexp +                    Schoolyears + ExpMilitary + Mortality +     PopPoverty + PopTotal + ExpEdu + ExpHealth, data=olympiadaten, nbest=2)     summary(a)     plot(a,scale="adjr2")               summary(lm(Gesamt~ExpHealth))      screenshot of the plot:   ![http://tinypic.com/r/2pq8agy/6](http://i.stack.imgur.com/GoRxx.png) The problem is now, that I want to fit the best model again "manually" and have a look at it, but the value of the adjusted R squared is not the same as in the regsubsets output? This is also the case for the other models, e.g. when I do the simplest model in the graphic:               summary(lm(Gesamt~ExpHealth))      The graphic says, it should have an adjusted R squared of about 0.14, but when I look at the output, I get a value of 0.06435. Here is the output of `summary(lm(Gesamt~ExpHealth))`:               Call:     lm(formula = Gesamt ~ ExpHealth)          Residuals:         Min      1Q  Median      3Q     Max      -18.686  -9.856  -4.496   1.434  81.980           Coefficients:                 Estimate Std. Error t value Pr(>|t|)       (Intercept)  -3.0681     6.1683  -0.497   0.6203       ExpHealth     1.9903     0.7805   2.550   0.0127 *     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 18.71 on 79 degrees of freedom       (4 observations deleted due to missingness)     Multiple R-squared: 0.07605,    Adjusted R-squared: 0.06435      F-statistic: 6.502 on 1 and 79 DF,  p-value: 0.01271       I don't know what I might have done wrong, any help would be appreciated. And last but not least, some more questions:   * What is the difference between selecting models by AIC and by the adj. R squared?    * Both measure the fit and recognize the number of variables, so isn't the best model chosen by AIC also the model with the highest adj. r squared?    * When I have 12 variables, this means, there are $2^12$ possibilities of models, right?    * So does the `regsubsets()` command calculate each model and show the two best (`nbest=2`) of each size?    * If so, do I really get the 'best' model?    * And when I do AIC using backwards selection (starting with the model which contains all variables), does this also end up with the same model that `regsubsets()` says is the best?


title: Poisson regression with both response and explanatory variables as counting

text: I’ve got the following variables: Response: number of quota units leased (in and out) by fishers. Explanatory: number of quota units own by fishers. I fitted a GLM (Poisson), but I’m not totally sure if it’s right, considering that the explanatory variable is count as well. I’ve found examples of Poisson regression just with categorical and continuous explanatory variables, but not with counting variables. So: Am I right using Poisson with my data? If not so, what alternative do I have? The residuals variances of my model are not homogeneous. I understand that Poisson regression allows face this problem, or should I pay attention to this issue and solve it (using weights, for example)? Any help would much appreciated, Cheers


title: How representative is Poisson distribution of the distribution of events in reality?

text: I've always wondered how good a 'fit' is the Poisson distribution to the events we observe in reality. Almost always I've seen it be used for modeling occurrence of events. (For example, arrival of cars in a parking garage or the number or messages sent/received by computers hosts on a network etc.) We usually model such events by the Poisson Distribution. Is the distribution just a good first approximation to how things happen in reality? If I observe the number of cars/day or messages/day in the above two examples and those that are output by 'picking from the distribution' how much do they differ? How good an approximation is Poisson? (Is it an approximation?) What is the 'magic' behind Poisson that it just gets it right (intuitively speaking :)?


title: Learning ML/NLP from electronics & communication engineering background

text: I study Electronics & Communication Engineering. I like Computer Sciences. Out of pure passion, I even studied few programming languages like Java, Python etc I studied probability, math, computer architecture, Java, C at college. I think Machine Learning, NLP etc stuff are fun for me. I never studied these subjects but I can say with the idea I have on them. So, can you please answer the below questions   1. I am even thinking to apply internships, can I apply in ML/NLP areas? What are the chances that I will get accepted?   2. Is a strong knowledge in ML/NLP is required to pursue UG internship in those areas? Is a strong desire to learn not enough?


title: May we categorize independent varibles per outcomes of dependent variables in logistic regression?

text: The question is rather theoretical than practical. Suppose, we want to evaluate one's chances of non-paying back a loan. We assume, that these chances depend on some set of covariates. One of the covariates is first two digits of ZIP code, let's call it zip2. We have a set of observations (loans) for some period of time, where for each loan we observe target variable y = 1, if loan is paid back, and y = 0 otherwise. Also, we have values of independent variable x = zip2. Let's say that # of observation is satisfactory, for example, 2000. Now, we are building a logit model in order to be able to predict the chances for the incoming loans to be (or not to be) paid back. Suppose, we run some software package, like R or Stata to estimate the model with our dataset. In the result we see, that our independent variable zip2 is not significant in the model. Then, we try to analyze our dataset more closely. We see that percentage of "bad" outcomes for zip2 = "11" and zip2 = "73" is virtually the same, say 15%. So, we decide to make one category out of these 2. We also examine the other values of the independent variables, and categorize them in a similar fashion, so, in total, we end up with 6-7 categories. Now, when we run estimation for the model with "categorized" variable zip2, we see that it becomes significant! So, my question is: **Is this correct to do such categorization?** I mean, we group independent variable by outcomes of dependent variable, and then we try to predict the outcome. It seems like a vicious circle, and I feel that it is not correct at all. This example is taken from a real situation, existing in a risk department of one company. If this is wrong, I just want to point out that such approach to model building is not acceptable and they're simply making wrong decisions. I already consulted a professor in statistics, and he said that this was wrong to group independent variables with respect to outcomes of dependent variable. But, I want to find some theoretical proofs, why such grouping is wrong. Is there any book or article, which explains inconsistency of such approach? Or, maybe it is consistent?


title: How to obtain predicted values from a gamm() using averaged coefficients (MuMIn)?

text: I want to extract the predicted values from a gamm() whose coefficients have been averaged using the package MuMIn, but I'm getting an error.               library(mgcv)     library(MuMIn)     set.seed(0)     dat <- gamSim(6, n = 100, scale = 5, dist = "normal")     fmm1 <- uGamm(y ~s(x0)+ s(x3) + s(x2), family = gaussian, data = dat, random = list(fac = ~1))          dredge(fmm1)->sel     confset.95p <- get.models (sel, cumsum(weight) <= .95)     avgmod<- model.avg(confset.95p) #           pred <- predict (avgmod, type="link", backtransform= TRUE)               Error in predict.averaging(avgmod, type = "link", backtransform = TRUE) :        'predict' for models '3', '7' and '4' caused errors     In addition: Warning messages:     1: In UseMethod("predict") :       no applicable method for 'predict' applied to an object of class "c('gamm', 'list')"     2: In UseMethod("predict") :       no applicable method for 'predict' applied to an object of class "c('gamm', 'list')"     3: In UseMethod("predict") :       no applicable method for 'predict' applied to an object of class "c('gamm', 'list')"      From the error message it appears I am doing more than one thing wrong. Can anyone help understand this? Thanks


title: Interpreting Logit Interaction Term Coefficients (continuous * categorical)

text: I have the following output from a logistic regression model.               Coefficients:                                  Estimate Std. Error z value Pr(>|z|)         (Intercept)                -2.6023448  0.0694696 -37.460  < 2e-16 ***     our_bid                     0.0039520  0.0007646   5.169 2.35e-07 ***     our_bid:zipcode10000:14849  0.0019334  0.0009006   2.147 0.031807 *       our_bid:zipcode14850:19699  0.0022905  0.0009514   2.407 0.016064 *       our_bid:zipcode19700:29999 -0.0009483  0.0008583  -1.105 0.269231         our_bid:zipcode30000:31999 -0.0016309  0.0011028  -1.479 0.139161         our_bid:zipcode32000:34999  0.0016241  0.0007856   2.067 0.038688 *       our_bid:zipcode35000:42999  0.0023549  0.0008541   2.757 0.005831 **      our_bid:zipcode43000:49999  0.0007096  0.0008104   0.876 0.381286         our_bid:zipcode50000:59999  0.0006533  0.0009269   0.705 0.480942         our_bid:zipcode60000:69999  0.0030564  0.0008169   3.742 0.000183 ***     our_bid:zipcode7000:9999   -0.0027419  0.0012699  -2.159 0.030847 *       our_bid:zipcode70000:79999  0.0013243  0.0007809   1.696 0.089921 .       our_bid:zipcode80000:89999  0.0038726  0.0008006   4.837 1.32e-06 ***     our_bid:zipcode90000:96999  0.0038746  0.0007817   4.957 7.18e-07 ***     our_bid:zipcode97000:99820  0.0009085  0.0010044   0.905 0.365726         ---      I am using these coefficients to draw the predicted probabilities such that. $$\text{Prob} = \frac{1}{1 + e^{-z}}$$ where $$z = B_0 + B_1X_1 + \dots + B_nX_n.$$ I realize that interpreting these interaction terms can be challenging. However, I generate the main regression equation and use that to formulate the probability curve. However, I'm not sure how to make sense of any of the "our_bid:zipcode" variables? What about if my model output was: (instead saving zipcode as a factor, I make it a continuous variable)               Coefficients:                                  Estimate Std. Error z value Pr(>|z|)         (Intercept)                -2.6023448  0.0694696 -37.460  < 2e-16 ***     our_bid                     0.0039520  0.0007646   5.169 2.35e-07 ***     our_bid:zipcode             0.0019334  0.0009006   2.147 0.031807 *        Would interpretation being easier with this approach? Keeping with the log- odds, how can I make sense of the log-odds effect that this model expresses for the interaction term?


title: Finding P Value of Two sample hypotesis test for Proportion

text: I am attempting to find the P-value of a hypothesis test, which compares the success rate of two procedures. A test with 0.05 level of significance is conducted to test the success rate of two procedures. **Procedure A** Sample Size: 200, Number of Success: 125 **Procedure B** Sample Size: 100, Number of Success: 55 **My Attempt** > ![enter image description here](http://i.stack.imgur.com/KecSK.jpg) The cooprrect p-val in the ans key is 0.2113. Appreciate some guidiance on where Ive gone wrong


title: Paired comparison for big data

text: I want to measure the subject's response to 30 stimulus. I'm thinking of paired comparison, but this experiment will be impossible to perform (need C(30,2)=435 comparison steps). The reason of using paired comparison is its accuracy (as opposed to simply ranking) What are my other alternative options? Can I sample the comparison space, and interpolate it?


title: Strategy for editing comma separated value (CSV) files

text: When I work on data analysis projects I often store data in comma or tab- delimited (CSV, TSV) data files. While data often belongs in a dedicated database management system. For many of my applications, this would be overdoing things. I can edit CSV and TSV files in Excel (or presumably another Spreadsheet program). This has benefits:   * spreadsheets make it easy to enter data There are also several problems:   * Working with CSV and TSV files leads to a wide range of warning messages about various features being lost and how only the active sheet will be saved and so forth. Thus, it's annoying if you just want to open the file and make a little change.   * It does many "supposedly intelligent" conversions. For example, if you enter 12/3, it will think that you want to enter a date. **UPDATE:** I should have mentioned that the date example is just one of many examples; most problems seem to be related to inappropriate conversion. In particular, text fields that look like numbers or dates cause problems. Alternatively, I could work directly with the text file in a standard text editor. This ensures that what I enter is what is recorded. However it is a very awkward way to enter data (columns don't line up; it's difficult to enter data simply into multiple cells; etc.). Question   * What is a good strategy for working with CSV or TSV data files? i.e., what strategy makes it easy to enter and manipulate the data while also ensuring that what you enter is actually interpreted correctly?


title: Visualizing the distribution of something within a very large body of data

text: Say I have a database of around a million words, and I want to get an intuitive idea about exactly how a particular, quite infrequent, word is distributed throughout this data. My goal is to be able to see clearly whether this word tends to cluster together, or whether it is relatively evenly spaced. What would be some good methods for visualizing this? For instance, I have seen something that looks useful. It's basically a strip (long rectangle) in which each instance of something is represented by a very thin red vertical line. The problem is that I don't know what these are called, and therefore I can't figure out how to make something like this in R. Any help finding the right R function for that, or any other suggestions for good ways to visualize this sort of data, would be most appreciated.


title: Method to compare variable coefficient in two regression models

text: I am regressing two butterfly richness variables (summer and winter) against a set of environmental variables separately. (variables with continuous numbers) Environmental variables are identitcal in each model. In the summer model, the weight rank of coefficients is temp > prec > ndvi. The weight rank in winter is temp > ndvi > prec. As it is almost implausible to compare the coefficients directly, pls advise any advanced method other than regression to discriminate such coefficient rank between seasons, such as canonical correlation analysis (unsure if it is suitable here) The spatial info here richness and environmental variables comprising 2000 grid (continuous distribution) spanning from 100 E to 130 E longitude, 18 to 25 N latitude.


title: Show that the cross-product term $\hat y'r$ from $ y'y$ is 0

text: Show that the cross-product term \hat $y'r$ from $y'y$ is 0. I see that $y'y = ( \hat y + r)'( \hat y + r) = \hat y' \hat y + \hat y'r + r' \hat y + r'r$ But not sure what identity can be used here to break up $\hat y'r$ such that it can be equal to 0. r = residuals y = response value $\hat y$ = predicted value The notes package written by a professor has this under an ANOVA section so the link you gave me, I'm not seeing the connection, but then again I also don't get the course notes and question so maybe I'm wrong here :|


title: How to find summary statistics for all unique combinations of factors in a data.frame in R?

text: I want to calculate a summary of a variable in a data.frame for each unique combination of factors in the data.frame. Should I use plyr to do this? I am ok with using loops as opposed to apply() ; so just finding out each unique combination would be enough.


title: Survival Analysis tools in Python

text: I am wondering if there are any packages for python that is capable of performing survival analysis. I have been using the survival package in R but would like to port my work to python.


title: Comparing proportions across more than 2 groups

text: I'm a clinical researcher and I am comparing children who had fluid in the abdomen (ascites due to liver disease). I have 4 groups based on the kind of fluid analysis and in each group I have incidence of renal failure, bleed, altered sensorium, etc. How should I compare these groups? Should I do multiple chi square tests taking 2 groups at a time? My number in each group is above 30 (but the number is different in each), and I am using SPSS 20.


title: Compare independent numeric vectors (means) in R

text: I am working on R, and have multiple data.frames. They all look like this one:                head(Stim.Adjs)           LogFreq    Word   PhonCV    FreqDev  LengthDev     6316    2.673  enkele VC-CV-CV -0.1794029 -1.0608453     12638   2.415 laatste CVVC-CCV -0.2207496 -0.7113483     12639   2.415 laatste CVVC-CCV -0.2207496 -0.7113483     1633    2.248    bang      CVC -0.2475127 -1.7598394     3167    2.243   bezig  CVV-CVC -0.2483140 -1.4103424     25467   2.243 vroeger CCVV-CVC -0.2483140 -0.7113483      These are list of words extracted from a corpus, and the columns 4 and 5 represent measures of standard deviation of Frequency and Word length respectively. What I want to do is to have all lists with a similar frequency deviation and word length deviation (two separate controls). So far I had been doing this by comparing each dimension in pairs of list using t.tests, and then removing values from the lists so that the t.tests showed no significant difference in means. Now that I have more than two lists, I'm a bit confused as to how to test if the means are equal or not. At first I though of an ANOVA, but I don't have a predictor variable, just 3 independent means (for each of the two values I want to control). The question is, then, what alternative to the t.test could I use that allows me to compare the means of independent vectors. Help with the function for this in R is much appreciated. Thanks in advance.


title: Relative density plots in R

text: I have a dataset with two classes, A and B, and an independent variable x. Class B is a minority class with only 1% of the observations, which tend to occur with small values of x. I am interested in the relative densities for the classes, which I would like to visualize in R. In my imagination, I would use `plot(density(data[cl=="A"]$x)/density(data$x))`, but R doesn't allow division of densities. How should I do this?


title: How to use the SD of a normal sampling distribution to specify the gamma prior for the corresponding precision?

text: The gamma distribution is a commonly used prior distribution for the precision ($1/sd^2$) of a normal distribution in Bayesian hierarchical modeling. I want to use an informed prior for the variance and I ask an expert on the data I'm trying to model. An easy way of describing earlier knowledge about the variance of the data is using the standard deviation as it is on the same scale as the data. My expert tells me: "Well in general the data tend to have a SD of 100, and I believe the SD of the SD is at most 50". I believe it would be much harder for my expert directly give this statement using precision and the SD of the precision. Now that I have the expert's statement I want to incorporate it into my JAGS/BUGS model but in JAGS/BUGS I can't specify the gamma prior directly using this data as gamma is parametrized using shape and rate and the normal distribution is parameterized using precision. What I would want to do is to take the statements of the mean of the sampling distribution's SD ($\mu_\sigma$) and the SD of the sampling distribution's SD ($\sigma_\sigma$) and use these parameters to set up the corresponding gamma prior for the precision of the sampling distribution. How could I do that? **That is, given the expert's guestimate of the mean SD and the SD of the SD of the data, how do I specify the corresponding hyperparameters shape and rate for the gamma prior on the sampling distributions precision?** The following diagram shows the model and illustrates my question: ![enter image description here](http://i.stack.imgur.com/MukIn.png) The code below would be the corresponding BUGS/JAGS model, where the "functions" `f1` and `f2` are the functions I would like to know how to implement.               model {       for(i in 1:length(y)) {         y ~ dnorm(mu, tau)       }            tau ~ dgamma(shape, rate)       shape <- f1(mu_sigma, sigma_sigma)       rate <- f2(mu_sigma, sigma_sigma)       mu_sigma <- 100       sigma_sigma <- 50            mu ~ dnorm(M, P)       M <- 0       P <- 0.0001     }      The reason why I would want to specify the shape and the rate of the gamma using $\mu_\sigma$ and $\sigma_\sigma$ is because I would find it more intuitive to think in the scale of the SD of the sampling distribution rather than in the scale of the precision.


title: what does histogram mean in a deterministic signal?

text: I have a deterministic signal comming from an experiment: x1[n]=[2,3,2,4,1,4] as a **set** ={1,2,3,4} histogram={1/6,2/6,1/6,2/6} what does the histogram(number of occurrence of each value) mean in the case of a deterministic signal such as X1[n]? It's not the probability density function, since we don't define a probability density function for a deterministic signal. Thanks for any clarification


title: Is there any guide line to select a right coefficient of variation?

text: I faced an estimation problem of compute the mean of a random variable (rv) $X$. What I know is the coefficient of variation of this rv, namely its standard deviation divided by its mean: $$c_{v}(X) = {\sigma_{x}\over \mu_{x}}$$ Surely, I know that the more samples I have, the more accurate my estimation of $\mu$, but I want to use as small number of samples as possible (while keeping my estimation reasonably sound). If we look the sample mean $\bar{X}=\sum_{i=1}^N{X_i}/N$, we know that $$c_{v}(\bar{X}) = {1\over \sqrt{N}} c_{v}(X)$$ My plan is to choose the number of samples $N$ with respect to $c_{v}(\bar{X})$, but I am not sure what a good $c_{v}(X)$ should be. Is there any rule of thumb? PS: $c_{v}(X)>1$.


title: Statistical significance problem

text: I've got an algorithm which assigns a score to a list of individuals (nodes). I want to discover which nodes received a statistically significant score. So I "randomized" the algorithm in order to obtain 100 random rankings of the same individuals. Now I've got my 100 rankings (Original + 100 Random). How can I check which nodes obtained a statistically significant score comparing the original with the random ones? Which statistical test could fit with my situation? EDIT: What is I define my null hypothesis as: "the observed value (original score) of node X appeared by chance". Then if the score of my original node X is higher than the 99 random scores, can I say that I reject the null hypothesis with a p-value < 0.01?


title: Mathematical model building with dependent and independent variables

text: I have been working with data and building models on data. I have developed models in regression using cubic and power series. It works fine for variables with one dependent and one independent variable. What kind of technique is to be adopted if I have to find an equation or model for say, D depends on C, C changes for a set of B, which changes for different A. I can carry our bi-variate correlation to find significance. But, that is not something I am interested in. I have data, but wondering as to how to fit a model in those lines. I have tried building models with regression, line fitting and get a relationship between D and C, or D and B. Good fits with varying constants, but no relationship between the others. Is there a way to combine and build a model altogether? What is that kind of an analysis called? How to go about it? Any softwares that would help me do it? I have been using SPSS, Minitab and R.


title: What are the issues when performing an ANCOVA with a covariate that is influenced by experimental manipulation?

text: ### Context: I am doing a study on lying.   * The independent variable was manipulated such that participants were either told to **lie** or **tell the truth**.   * The dependent variable is related to the amount of detail provided and is extracted from the text provided by participants(this is a standard system of deception detection called reality monitoring).  ### Issue: I noticed that much of the previous research (e.g., Sporer [2000]) has included or suggested the use of word count as a covariate to make sure that any differences in the number of details do not simply reflect differences in the number of words. However, I was under the impression that covariates needed to be measured before the experimental manipulation and should not be related to the independent variable. Word count is measured after people have provided their truthful or deceptive account (And possibly related to truthfulness since people might provide more details by giving a longer account). I have read Miller and Chapman's paper on the misunderstanding of ANCOVA (FREE PDF). ### Questions   * Does including word count as a covariate looking at the effect of condition on details provided reduce statistical power?   * Is it meaningful to include a covariate that is measured after, and is assumed to be affected by, an experimental manipulation?   * How serious is violation of heterogeneity of covariance when including covariates?   * Rather than including word length as a covariate, often authors calculate number of details say per 100 words and then perform their analyses as normal. how is this different to including word length as a covariate? ### Additional References: For those interested, references include:   * Stromwall, L., & Granhag, A. (2005). Children's repeated lies and truths: effects on adults' judgements and reality monitoring scores. Psychiatry, Psychology and Law, 12(2), 345-356,    * Vrij, A., Mann, S., Fisher, R., Leal, S., Milne, R., & Bull, R. (2008b). Increasing Cognitive Load to Facilitate Lie Detection: The Benefit of Recalling an Event in Reverse Order. Law and Human Behavior, 32(3), 253-265.   * Sporer, S. L. (2004). Reality monitoring and the detection of deception. In P. A. Granhag & L. Stromwall (Eds.), The detection of deception in forensic contexts. Cambridge: Cambridge University Press.


title: Is it possible to graph this?

text: Is it possible to graph this? I tried doing it on a calculator but didn't have any success. If someone can and post the graph I would greatly appreciate it. $$\frac{1}{\bar{m}_{T}}\sqrt{\frac{1}{12} \sum_{h=T-11}^{T} (m_{h} - \bar{m}_{T})^2}.$$ If the equation is hard to read, it is just the standard deviation of 12 observations $(m_{h})$, $h=T-11, T-10, \ldots, T$, over the mean of the same 12 observations. The standard deviation is the population standard deviation and the mean $\bar{m}_{T}$ is simple average.


title: Grouping data in ranges in R by summing them

text: I have data like this :                     year   nb     1     1901  208     2     1902  200     3     1903  223     4     1904  215     5     1905  187     6     1906  214      And I want to specify levels, such that I can summarize the data this way :                     years   nb     1     1901-1910  2082     2     1911-1920  6200      I had a hard time doing this either with group, aggregate, or encode until then. I found a very ugly way of doing that, like that :               sum(DF$nb[DF$year> 1901 & DF$year <= 1910])      But I would like to know if there is a more elegant way to do it. Sorry if my question is too basic, Xavier


title: Linear regression for samples of varying sizes

text: I have a data set of code patches and the bugs they produce. I'm using ordinary least squares to find a line which predicts bugs based on some attributes about the patch, such as the department which produced the code, how many issues were found in design, etc. One thing is that the patches vary in size, and larger patches will (obviously) create more bugs. I'm trying to think of how I should handle this, and I came up with several options:   1. Use patch size (lines of code changed) as an independent variable. The problem with this is that I think I really want to predict defect "density". E.g. the Sales department might cause one bug every 15 lines of code, but the shipping people cause one every 30 lines of code. So I would need a variable like "department * lines of code changed", which I worry is more sensitive to overfitting.    2. Break it into samples of equal size. E.g. if a change to 10 lines of code caused 2 bugs, I pretend that this was actually 10 changes, each of which caused .2 bugs.   3. Divide through by the size and then weight by the size. i.e. minimize $\sum w_i(\frac{y_i}{w_i} - f(\vec{x_i}))^2$.  Are any of these preferable, or is there another option I should be using? I assume this must be a pretty standard problem to have, but I haven't been able to find anything online.


title: Interpreting as much data from a single statistical value

text: A song has a total of 1,000 plays, and on average only 50% of the song was played (say a song is 4 mins long and on average listeners only listen to half the song). What other data can we interpret from the numbers above? I'm trying to estimate how many of the 1,000 plays actually listened to the whole song, but it doesn't seem possible, does it? Thank you.


title: How to use a mathematical model for data analysis in R

text: I am looking to use a mathematical model developed by Firbank & Watkinson (1985) J. App. Ecol. 22:503-517 for the analysis of competition between plants grown in mixture. The model is as follows: $$W_{A}=W_{mA}\left(1 + a_{A}\left(N_{A}+\alpha N_{B}\right)\right)^{-b_{A}}$$ where $W_{A}$ is the mean yield per plant of species $A$ grown in the experiment, $W_{mA}$ is the mean yield of isolated plants of species $A$, $a_{A}$ is the surface area required to reach size $W_{mA}$, $N_{A}$ is the planting density of species $A$, $\alpha$ is the competition coefficient, and $-b_{A}$ is the 'resource use efficiency' parameter. The model is a regression model as I understand it. I have data for density of species $A$ and $B$ and ($N_{A}$ and $N_{B}$) as well as the response variable $W_{A}$. I am unsure how I can use R to estimate the remaining values, most important of which is the competition coefficient, $\alpha$. If there is any more information that I need to provide please let me know.


title: ACF of IMA(1,1) in R

text: I generate the IMA(1,1) process in R to see ACF for different values of Theta. The graphs show that ACF is always positive for both positive and negative values of Theta, is it correct or I'm wrong somewhere?               theta <- 0.9     d <- ts(matrix(0,ncol=1,nrow=100))     e <- ts(rnorm(100,0,1))     d[1] <- e[1]     for(i in 2:100)       d[i] <- d[i-1]-theta*(e[i-1]+e[i])     acf(d,10)      Theta=0.9 ![enter image description here](http://i.stack.imgur.com/pPyE6.png) Theta=-0.9 ![enter image description here](http://i.stack.imgur.com/vnlvp.png)


title: Marketing Mix Modelling

text: I am looking for some text book type of explanation how the Marketing Mix Modelling works, and what Statistical tools are used in doing the same. Hands on examples with explanation would really be good. Does anyone come up with any text book on this subject? It would be really helpful if people can share their ideas on the same. Thanks for your inputs.


title: Hyper-prior for negative binomial in hierarchical model using JAGS/BUGS

text: Below I'm using a negative binomial because it is more flexible than a simple poisson model. The data are counts $y$ of events for 16 individuals $x$. There are 14 counts (i.e. counting periods) for each individual. The likelihood function is `dnegbin`($p_x , r_x)$ with a separate $p$ and $r$ for each individual $x$ drawn from a gamma distribution. (Actually I have re- parameterized this bit according to Kruschke.) I follow Kruschke's choice of gamma priors for $r$ and $m$ here. Now my question is how to pick the hyper-priors for the parameters of the two gamma distributions from which the parameters of the negative binomial come:               r[j] ~ dgamma(r.a,r.b)     m[j] ~ dgamma(m.a,m.b)      I've experimented with gamma distributions `dgamma(a,b)` and setting $a$ and $b$ both to 0.5 or 1.0 seems to yield reasonable results.   1. Is the choice of gamma priors for both parameters of a gamma distribution reasonable, i.e. should both `r.a` and `r.b` come from a gamma distribution?   2. What are the parameters of a minimally informative gamma prior in this case, e.g. `dgamma(0.5,0.5)`? . # The JAGS/BUGS model:               modelString = "     model {         for( i in 1:length(y) ) {           y[i] ~ dnegbin( p[x[i]] , r[x[i]] )         }         for( j in 1:Nj ) {             p[j] <- r[j]/(r[j]+m[j])             r[j] ~ dgamma(r.a,r.b)             m[j] ~ dgamma(m.a,m.b)             v[j] <- r[j]*(1-p[j])/(p[j]*p[j])         }         r.a ~ dgamma(a, b)         r.b ~ dgamma(a, b)         m.a ~ dgamma(a, b)         m.b ~ dgamma(a, b)         a <- 1.0         b <- 1.0     }     "      . # Running the model For experimentation I provide some data:               d <-     list(x = rep(1:16, each=14),     y = c(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 1, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 4, 0,     2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 2, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,     3, 0, 0, 0, 0, 0, 4, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 2, 0,     0, 3, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0,     1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 0, 1,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), Nj = 16)      Run JAGS:               require(rjags)     parameters <- c("m","r","p","v","r.a","r.b","m.a","m.b")     jagsModel <- jags.model(textConnection(modelString), data=d, n.chains=3, n.adapt=1e3)     update(jagsModel, n.iter=1e3)     codaSamples <- coda.samples(jagsModel, variable.names=parameters, n.iter=3000, thin=1)     m <- as.matrix(codaSamples)      Extract parameters (random effects) for each individual:               out <- data.frame()     for (p in c('m','r','p','v')) {         for (i in 1:16) {             out[i,'x'] <- i             c <- sprintf('%s[%i]',p,i)             #out[i,p] <- mean(estimates[c])        # mean             dens <- density(m[,c])             out[i,p] <- dens$x[which.max(dens$y)]  # mode         }     }          # calculate frequencies for each individual 'x' and #occurrences      t <- data.frame(table(d$x,cut(d$y, breaks=0:6, right=F)))     colnames(t) <- c('x','breaks','freq')          # add fits to frequency table     t <- merge(t,out)     t <- within(t, freq.fit <- dnbinom(as.numeric(breaks)-1, size=r, mu=m)*14)      Plot results:               library(ggplot2)     ggplot(t,aes(breaks,freq,group=x)) + geom_point() + geom_line(aes(y=freq.fit)) + facet_wrap(~x) + theme(panel.grid=element_blank(), panel.background=element_blank())      For anyone who is struggling to relate the arguments of the negative binomials from JAGS and R:               JAGS,  R's dnbinom        r,         size        p,         prob        m,           mu


title: SVM with unequal group sizes in training data

text: I am trying to build an SVM from training data where one group is represented more than the other. However, the groups will be equally represented in the eventual test data. Therefore, I'd like to use the `class.weights` parameter of the `e1071` R package interface to `libsvm` to balance the influence of the two groups in the training data. Since I was unsure exactly how these weights should be specified, I set up a little test:   1. Generate some null data (random features; 2:1 ratio between group labels)   2. Fit an svm with the `class.weights` parameter set.   3. Predict a bunch of new null datasets and look at the class proportions.   4. Replicate the whole process many times for different null training sets. Here is the R code I'm using:               nullSVM <- function(n.var, n.obs) {         # Simulate null training data         vars   = matrix(rnorm(n.var*n.obs), nrow=n.obs)         labels = rep(c('a', 'a', 'b'), length.out=n.obs)         data   = data.frame(group=labels, vars)              # Fit SVM         fit = svm(group ~ ., data=data, class.weights=c(a=0.5, b=1))              # Calculate the average fraction of 'a' we would predict from null test data         mean(replicate(50, table(predict(fit, data.frame(matrix(rnorm(n.var*n.obs), nrow=n.obs))))[1])) / n.obs     }          library(e1071)     set.seed(12345)     mean(replicate(50, nullSVM(50, 300)))      From this whole thing I was expecting an output ~ 0.5, however, that's not what I got:               > mean(replicate(50, nullSVM(50, 300)))     [1] 0.6429987      The `class.weights` paramter is working, _sort of_ , as the lower I weight `a`, the lower it is represented in this simulation (and if I omit `class.weights` it returns close to 1)...but I do not understand why simply using weights of 1:2 (for training data that is 2:1) does not get me all the way down to 50%. If I'm misunderstanding SVMs, can someone explain this point? (or send some refs?) If I'm doing it wrong, can someone tell me the correct way to use the `class.weights` parameter? Could it possibly be a bug? (I think not, since I understand this software and the underlying libsvm to be quite mature)


title: Probability of circles intersecting

text: Suppose that N circles are dropped randomly onto a square in 2-space in such a way that the center of each circle lands somewhere in the square – even if some of the circle extends beyond the boundary of the square. The side of the square is of length S, where S is an integer. Each circle has a real number radius R, where R is a uniformly distributed random value between 0 and some integer M. (Each circle has its own radius of size between 0 and M units). What is the expected number of circles that will not intersect with any other circles? That was a question from a programming competition 3 weeks ago (http://www.cs.rit.edu/~icpc/questions/2011/Oswego_2011.pdf). My team couldn't solve it, but I really want to know how to. It even inspired me to take a probability course next quarter because I want to know more about that subject.


title: How to open ws. files in datasets?

text: I have to use the datasets on http://www.bristol.ac.uk/cmm/learning/support/datasets/. And I found every data set on this website has a ws. file. I try to download software which can open it and I can get the data. But I do not find a free one. It is very disappointed. Could you recommend one for me?


title: Logistic regression with block correction for ordinal response

text: I have a categorical response (ordinal, lets call it A) collected from a complete randomized block design experiment where the explanatory variable is treatment. How do I account for the between block variation? I have found models that do this for dichotomous response variables but my response has >2 levels, treatment is also categorical and has 10 levels and there are 8 separate blocks. Could I hypothetically run a logistic regression and add the separate block factors to account for the variation? The treatment variable has 10 levels but I could separate them so that I compare them 2 at a time, control vs some treatment. Would it be reasonable to run a GLMM model using treatment as a response? So that level changes in the categorical variable, A, affect the odds of having received treatment or is there some fallacy in this?


title: Is it possible to fit a multilevel Weibull/generalised gamma survival model?

text: I am imagining a situation in which I have individual patient survival data from a number of clinical trials of a particular drug. The drug can be used for a wide variety of cancers which have significantly heterogeneous survival distributions, but suppose that a Weibull or generalised gamma distribution would be an appropriate fit in each case. The data will include right-censored patients, so the data could look like:               PatientID  StudyID  LastObservation  Died  Treatment  Male     ----------------------------------------------------------     1          1        3.22             1     0          1     2          1        5.00             0     0          0     3          1        6.49             1     1          0     4          1        0.77             1     0          1     5          1        2.50             0     1          1     ...      In this example data `Died = 0` means the patient is right-censored and `Male` is included as a covariate (in reality there would be many more covariates, e.g., age at randomisation, cancer site, baseline haemoglobin level). I would like to fit a model to this, e.g., $f\left(t_{ij}\right) = \lambda_{ij}\gamma_{ij}t_{ij}^{\gamma_{ij}-1}\exp\left(-\lambda_{ij}t_{ij}^{\gamma_{ij}}\right)$ where $\ln{\lambda_{ij}} = \pmb{\alpha_j}\cdot \pmb{u_{ij}}$, and $\ln{\gamma_{ij}} = \pmb{\beta_j}\cdot \pmb{u_{ij}}$, with $\pmb{\alpha_j}$ and $\pmb{\beta_j}$ drawn from multivariate normal at the study level $j$ and $\pmb{u_{ij}}$ is the vector of covariates (including `Treatment`) and a constant term 1. This would accommodate random intercepts and slopes as necessary, e.g., the impact of treatment on survival can vary between studies. Is it possible to fit such a model? Probably doesn't need to have a full covariance matrix (block diagonal may suffice). I would be happy to use generalised linear (multilevel) model framework if appropriate. My software options are R, STATA, MLwiN and WinBUGS.


title: Interpretation of a WEKA result buffer - confusion matrix and performance

text: I want to know how to get several performance measurements of a generated WEKA model. Note that I am predicting a two-class variable, `Alive` or `Dead` and I use the developer version `3.7.10` in my own Java project.   * Confusion matrix I know this is a simple question but am I correct that the one below is a direct equivalent of a traditional confusion matrix?              === Confusion Matrix ===               a     b   <-- classified as      13735   216 |     a = Alive        392   657 |     b = Dead      As in:              TP | FN     FP | TN        * _Sensitivity_ and _specificity_ Given a part of the WEKA result buffer below, below contains the ROC, specificity (or `Recall`) and sensitivity (or `Precision`) measurements but depending on a given class.              === Detailed Accuracy By Class ===                           TP Rate  FP Rate  Precision  Recall   F-Measure  MCC      ROC Area  PRC Area  Class                      0.985    0.374    0.972      0.985    0.978      0.665    0.974     0.998     Alive                      0.626    0.015    0.753      0.626    0.684      0.665    0.974     0.775     Dead     Weighted Avg.    0.959    0.349    0.957      0.959    0.958      0.665    0.974     0.982           Poor Felix asked the questions " _WEKA Specificity and Sensitivity for global Rule Model_ " in (a), (b) and (c) similar to the one I have right now but still didn't get the exact answer. Now, is it right to interpret that the overall ROC, specificity and sensitivity of the whole model are the ones in the `Weighted Avg.` row?


title: IVEware imputation results

text: I have a large survey dataset (~150k individuals) on which I used IVEware to impute missing data. Most of the variables missing data are missing very small percentages of data (<5%). However, many time these missing values are on different cases, so complete case analysis would result in dropping a far greater percentage than 5%, hence the motivation for imputation. I'm considering about 40 variables. On about 4 of these , IVEware imputed a single value, so there is no between imputation variance. These question is what to do about it.               1) Should I go ahead an use the imputed data and note that these variables did not have      between imputation variance?     2) Is this a complete failure leaving my data useless?  What then?  I ran 10 imputations with      20 iterations (at great expense of time).  Should I run more imputations and fewer iterations      to try to obtain the between imputation variance?  My worry here is that it will produce a      large file that will take a lot of time to process.


title: PCA, non positive definitive & potential cures?

text: I am trying to use PCA to remove the correlation in my dataset My correlation matrix is being reported as non positive definitive As I understand it a NPD report means that my transposed matrix does not contain all positive eigenvalues values. (i know what eigenvalues represent, but not how they are calculated) This maybe due to very strong near/perfect colinearity between 2 or more variables...this is quite possibly the case, I have removed all very high correlations but understand colinearity may still be present, this is another highly correlated dataset. Using row wise method does remove the NPD warning. I also understand this error may be due to more variables than cases. I have ~10x more cases than variables. I also gather that depending on the estimation methods cases can be dropped dependent on missing values. I do have some missing values but do not think there are enough to make cases less than variables Despite the warning I proceeded to perform PCA on both pairwise an row wise method results. With both methods the unrotated factors when put into a model do remove all high VIF issues I had before PCA. But as you would expect the factors themselves are not very interpretable. On rotation (varimax) the NPD pairwise method results in rotated factors that are very interpretable but colinearity is not removed. is this related to NPD warning? Again, after rotation (varimax) with the row wise method the results were no longer colinear at all. BUT the rotated factors were very much less interpretable as real world variables. I measured colinearity in the factors by putting all into a linear regression model and looking at VIF scores What options do i have here?


title: recommenderlab, why does the similarity function return only 0's and 1's?

text: In a matrix with user ratings I want to calculate similarities between the first ten users and the rest of the users. I use the method "jaccard" for the similarity calculation in recommenderlab. My matrix looks like this:                  users v1 v2 v3 v4 v5     1      1  0  1  0  0  1     2      2  0  1  0  0  1     3      3  0  1  0  0  1     4      4  0  1  0  1  1     5      5  0  1  0  1  1     6      6  0  1  0  1  0     7      7  0  1  1  1  1     8      8  1  1  1  1  1     9      9  0  1  0  1  1     10    10  0  1  0  1  1     11    11  0  1  1  0  1     12    12  0  0  0  1  1     13    13  0  1  0  1  1     14    14  0  1  1  1  1     15    15  0  1  0  0  0     16    16  0  1  1  0  1     17    17  0  1  0  0  1     18    18  0  1  0  1  1     19    19  1  1  0  1  1     20    20  0  1  0  1  1      First I make two datafiles               datafile1 <- datafile[1:10,]     datafile2 <- datafile[2:20,]      Now I try to calculatie similarities with the similarity function               library(recommenderlab)     datafile1 <- as(datafile1, "realRatingMatrix")     datafile2 <- as(datafile2, "realRatingMatrix")     similarity(datafile2,datafile1, method = "jaccard", which="users")      This is the result:                  1 10 2 3 4 5 6 7 8 9     11 0  1 0 0 1 0 0 0 0 1     12 0  1 0 0 1 0 0 0 0 1     13 1  0 0 1 0 0 0 0 0 0     14 0  0 1 0 0 1 0 1 0 0     15 0  0 0 0 0 0 1 0 0 0     16 0  1 0 0 1 0 0 0 0 1     17 0  1 0 0 1 0 0 0 0 1     18 0  1 0 0 1 0 0 0 0 1     19 0  0 0 0 0 0 0 0 1 0     20 1  0 0 1 0 0 0 0 0 0      I think this is not what I expected because I only see 0's and 1's.   Am I doing something wrong here? Many thanks in advance!


title: Normalizing variables '0_1' looks better than z-scores in R

text: Here is a toy-example, where '0_1' looks better than 'norm'. Using 'norm'(z-scores) has a huge training error and almost constant prediction.               library(RSNNS)     x <- runif(1020, 1, 5000)     y <- sqrt(x)     #0_1 type     x.01 <- normalizeData(x, type='0_1')     y.01 <- normalizeData(y, type='0_1')     model.01 <- mlp(x.01[1:1000], y.01[1:1000], size = 1)     py.01 <- predict(model.01, matrix(x.01[1001:1020], ncol=1))     py.01 <- denormalizeData(py.01, getNormParameters(y.01))     #norm type     x.n <- normalizeData(x, type='norm')     y.n <- normalizeData(y, type='norm')     model.n <- mlp(x.n[1:1000], y.n[1:1000], size = 1)     py.n <- predict(model.n, matrix(x.n[1001:1020], ncol=1))     py.n <- denormalizeData(py.n, getNormParameters(y.n))     #print     res <- cbind(y[1001:1020], py.01, py.n)     colnames(res) <- c('target', '0_1', 'norm')     print(res)         target      0_1     norm      [1,]  7.847747 18.42934 47.98046      [2,] 20.128710 21.63094 47.98193      [3,] 51.883655 54.31361 50.05302      [4,] 39.988633 38.61841 48.02621      [5,] 37.763832 35.78068 48.00804      [6,] 36.182139 33.88958 48.00035      [7,] 49.259397 51.15189 48.77056      [8,]  9.989678 18.74624 47.98058      [9,] 37.053090 34.91668 48.00420     [10,] 67.759140 65.00623 64.26272     [11,] 40.238815 38.94849 48.02908     [12,] 55.269081 57.80950 54.15643     [13,] 68.659288 65.26538 64.33199     [14,] 67.240699 64.84575 64.21361     [15,] 32.216679 29.69490 47.99008     [16,] 37.735459 35.74576 48.00787     [17,] 43.830574 43.84002 48.10858     [18,]  8.959444 18.58314 47.98052     [19,] 22.715183 22.83997 47.98266     [20,] 58.814483 60.69753 59.91301      Why doesn't using z-scores work here? I'm use RSNNS package


title: Axes on PCA-Biplot

text: When you construct a biplot for your PCA-analysis, you have comp.1 on the x-as and comp.2 on the y-as. But what are the other two axes to the right and the top of the screen?


title: ARIMA Models in R, implementing in C

text: I need to take the output parameters from an ARIMA model fitted in R from the following set (1,0,1), (0,1,0), (1,1,0), (0,1,1), (1,1,1) of models and implement the prediction function in C. I DO NOT HAVE THE OPTION of calling predict or any other R package for that step. Obviously, I can eventually track down all the source code in predict and figure it out. But I was hoping there is somewhere that will walk me through a simple example of how to map the various output parameters of Arima() with X, Y, a, b, E, t (no upper and lower case thetas and B^t's) since every paper loves to include those already. I think this request is slightly duplicative except in previous versions the question was retired without an answer or a link. UPDATE: So, first, I HIGHLY second all recommendations for Forecasting: Principles and Practice by Hyndman&Athanasopoulos. I think what I've been missing is that "d" isn't a model parameter -- it changes what is being modeled. So while I'm not all the way to where I want to be, I'm starting to be able to write predictive equations based on R output. I will update with my eventual findings if nobody else posts something better.


title: Regarding the formula of using $\text{P}(Y|X)$ to compute $\text{E}[X]$

text: When reading a presentation on "expectation propagation," I found a strange formula for computing $\text{E}[X]$ from a conditional probability: $$\text{E}[X] = \frac{\int x P(y_i|x) dx}{\int P(y_i|x) dx}$$ Can someone please explain how this formula was derived? The source is slide 24 in the presentation linked above.


title: How to report result from two-way repeated measures ANOVA?

text: I have conducted an empirical study as part of my master thesis. I have developed a software tool where cards can be sorted (affinity diagram) and in the study each test subject did a number of sorting tasks, to examine how the digital search feature improved completion time in sorting tasks. I am wondering if a Two-way, repeated measures ANOVA is appropriate to use when analyzing the data. I had two independent variables:   1. dataset size (10, 30 or 50 cards)   2. digital search feature ("on" or "off"). Which means I had 6 task in total, and the subject had to find 10 cards in each task, which was more or less difficult based on the dataset size and if digital search were enabled or disabled. I have two hypotheses that I want to validate:   1. Users will be able to achieve significant performance increases in card sorting tasks when using digital search.   2. Digital search are expected to play a small role when organizing a small set of cards. Digital search will therefore primarily outperform non digital search when organizing large sets of cards. It was a within-group study, with 14 participants (i.e. all subjects did all 6 tasks). The output is as follows (completion time in seconds):               |---------------No search--------------|----------------Search----------------|     |  10 cards  |  30 cards  |  50 cards  |  10 cards  |  30 cards  |  50 cards  |     |-----------------------------------------------------------------------------|     |     34     |    103     |     171    |     22     |     22     |     26     |     |     24     |     41     |      78    |     20     |     28     |     28     |     |     37     |     60     |     141    |     19     |     23     |     24     |     |     33     |     69     |     122    |     30     |     26     |     24     |     |     26     |     52     |     227    |     22     |     29     |     38     |     |     33     |     57     |     100    |     26     |     35     |     34     |     |     33     |     87     |     148    |     25     |     25     |     26     |     |     30     |     86     |     113    |     25     |     26     |     37     |     |     30     |    127     |     156    |     26     |     27     |     28     |     |     23     |     42     |     130    |     17     |     16     |     17     |     |     22     |     75     |     112    |     22     |     34     |     36     |     |     36     |     95     |     208    |     24     |     23     |     36     |     |     30     |     89     |     105    |     23     |     21     |     25     |     |     55     |    118     |     216    |     37     |     49     |     53     |     |-----------------------------------------------------------------------------|      I have run a Two-way, repeated measures ANOVA test in matlab with this matlab implementation which gave the following output:               'Source'                           'SS'          'df'  'MS'          'F'         'p'              'Dataset Size'                     [4.2174e+04]  [ 2]  [2.1087e+04]  [ 62.2973]  [1.2110e-10]     'Task type'                        [1.0926e+04]  [ 1]  [1.0926e+04]  [ 31.0805]  [8.9967e-05]     'Dataset Size x Task type'         [1.0708e+05]  [ 2]  [5.3539e+04]  [109.0535]  [2.2704e-13]     'Dataset Size x Subj'              [8.8008e+03]  [26]  [  338.4936]          []            []     'Task type x Subj'                 [4.5699e+03]  [13]  [  351.5311]          []            []     'Dataset Size x Task type x Subj'  [1.2765e+04]  [26]  [  490.9460]          []            []      My statistical knowledge is quite limited, so I am not sure if this significance test is relevant, or how to report the output. I guess the following can be used, when arguing that both hypotheses are true?   * F (2, 26) = 62.297, p < .001.   * F (1, 13) = 31.080, p < .001.   * F (2, 26) = 109.053, p < .001. * * * **Bonus questions:** Is a paired-sample t test appropriate to test if it took significant longer time to find 10 cards among 50 cards vs 10 cards when search were enabled (i.e. compare 4th and 6th column). I get a _significant_ difference (t(13)=3.76, p<0.01) if I use a t test.


title: Degrees of freedom due to treatment

text: let $a$ be the number of treatments. In Randomized Complete Block Design(RCBD), degrees of freedom due to treatment is $(a-1) $ . let $k$ be the number of treatments exactly in each block in a Balanced Incomplete Block Design(BIBD). In BIBD, Adjusted Mean Square due to treatment is given by adjusted sum of square due to treatment divide by degrees of freedom $(a-1) $. But my question is why do we divide it by $ (a-1) $? .Since it is adjusted , why don’t we divide it by $(k-1) $?


title: What is the null model for a likelihood ratio test of a within-subjects factor?

text: Tissue samples were taken from 4 differention locations and repeatedly measured. This was done identically for 3 animals. The research question was: Are there differences in measurement between the locations? => Repeated-measures design, with location as within-subjects factor               m1 <- lmer(meas ~ location + (1 + location | animal))        I would like to test for a main effect of location by a likelihood ratio test via anova(model, nullmodel) (I am aware of the differing opinions regarding this). But I am unclear, which, or if any, of the following two is the correct null model:               null.1 <- lmer(meas ~ 1 + (1 | animal), REML=FALSE)     null.2 <- lmer(meas ~ 1 + (1 + location | animal), REML=FALSE)      I am inclinced to think it’s null.1, but I am not entirely sure. Formulated more generally: When specifying a maximal random effects structure for a within-subjects factor, what is the null model for that factor for a maximum-likelihood ratio test?


title: Feature selection and training on the same sample

text: Is feature selection and training on the same sample a bad idea? I want to emphasize that I am not going to use test set for feature selection. If I use the whole train set for feature selection and then for training the model then the training sample is larger. Any drawbacks of such approach?


title: mixed ancova or difference scores?

text: I am puzzled with how I should analyze the following data: The dependent variable is respiratory sinus arrhythmia, a measure of parasympathetic activity (RSA). I have a $5\times5$ full factorial design: There are 5 levels of imposed breathing and 5 different moments in time. There are two groups (males and females) and a sample size of $n=20$. Additionally data from respiratory frequency (Fresp) in all conditions. Respiratory frequency has a direct effect on RSA, so I want to include this as a covariate. I.e. each data point of RSA contains information on Fresp. Apparently, `SPSS` does not allow that since the covariate is implied to be something that is related to the group differences. The possible approaches I have considered:   1. Standardize the RSA and the Fresp and subtract the standardized score from each other (so that I obtain variability in RSA 'free' of the variability of Fresp) and the perform a mixed anova on these data.   2. Perform a regression of Fresp on RSA within each individual and subject the residuals to a mixed anova.   3. Perform a multilevel regression with dummy variables and Fresp as independent variable. Are any of the approaches recommended or is there are different approach that someone would propose?


title: Test to distinguish periodic from almost periodic data

text: Suppose I have some function $f$ fulfilling some reasonable conditions like continuity. I know the exact values of $f$ (because the data comes from a simulation) at some equidistant sampling points $t_i=t_0 + iΔt$, which I can assume to be sufficiently fine to capture all relevant aspects of $f$, e.g., all local extrema, turning points, etc. I am looking for a test that tells me whether my data complies with $f$ being exactly periodic (for a reasonable period length). From another point of view, I have data ${x_0, …, x_n}$ and am looking for a test that answers the question whether a periodic function $f$ (fulfilling some basic conditions) exists such that $f(t_i)=x_i ∀ i$. The important point is that $f$ is at least very close to periodicity (it could be for example $f(t) := \sin(g(t)·t)$ or $f(t) := g(t)·\sin(t)$ with $g'(t) ≪ g(t_0)/Δt$) to the extent that changing one data point by a small amount may suffice to make the data comply with $f$ being exactly periodic. Thus standard tools for frequency analysis such as the Fourier transform or analysing zero crossings will not help much. Note that the test I am looking for will likely not be probabilistic. I have some ideas how to design such a test myself but want to avoid reinventing the wheel. So I am looking for an existing test.


title: How to test proportions of a Categorical response variable for a Repeated measures design with Unequal Sample Size?

text: I have a question about the analysis of a unequal sample size repeated measures data with categorical response variable. This experiment looks at 5 young and aged persons and for each of the individual we have a few data-points measuring their responses to a stimulus. The responses are categorical with four factors say A,B,C,D. So, for young person 1 we have $n^1_1$ observations, say (A,B,B,B,A,C,C, ...) and for old person 1 we have $n^2_1$ observations, say (C,D,D,D,C, ...,D). (Please note that $n_i^1$ and $n^2_i$'s might be very different). My goal is to compare the proportions of (A,B,C,D)'s among young people and older people. I am worried that an usual $\chi^2$ test might not be appropriate as the within-subject responses aren't independent. Can anyone suggest what would be the best way to go about testing the equality of proportions?


title: How to shape data into something useful for machine learning algorithms

text: I'm very new to this--I took a class last quarter and got hooked--now I'm in this competition and I don't know where to start (kind of). My problem is that I got this data that looks like a bunch of relational database tables and I don't know how to shape them into something I can try to feed into some kind of model. For example, I have users, users that are friends with other users, events, and users who are attending events. I know how to merge users and events together, convert it to floating point or integer values, and throw it into the model. Its the friends part that's has me confused. How can I bring the friends aspect into the equation?


title: Two seasonal periods in ARIMA using R

text: I'm currently using R to predict a time series with these instructions:               X <- ts(datas, frequency=24)     X.arima <- Arima(X, order=c(2,1,0), seasonal=c(1,1,1))     pred <- predict(X.arima, n.ahead=24)     plot.ts(pred$pred)      As you can see I've data each hour, and I chose the seasonal period of 24 (one day). I would like to improve my forecasting using an additional seasonal period in order to include the seasonal component of the week (seasonal length of 7*24=168 data) Is there any method for this? How do you do it? **UPDATE:** I've read this (your) blog page, maybe can I use the external regressors to simulate a second seasonal period?


title: What is the difference between an estimator and a statistic?

text: I learned that a statistic is an attribute you can obtain from samples.Taking many samples of same size, calculating this attribute for all of them and plotting the pdf, we get the distribution of the corresponding attribute or the distribution of the corresponding statistics. I also heard that statistics are made to be estimators, how do these two concepts differ?


title: Determine likelihood that a list of N items contains no negative results after examining X items

text: I have a large (N) list of results from a process that is either positive or negative. After examining X of them and receiving only positive results, how certain can I be that the list only contains positive results (with probability P)? Context: I'm processing the results (between several thousand and several billion) of a process and am given a list of values that take a few seconds each to evaluate as passing/failing. I only need to be reasonably certain that everything passed to pass the batch for most jobs. For some jobs, 100% certainty is necessary, and they'll all need to be checked. For the others, how many do I need to check to be P% certain that there are no fails?


title: How to determine influence of two time series with feedback?

text: Say you have two time series, A & B. Each have a mutual effect on the other. To give a real-world example, say that time series A measures an artist's CD sales per month, and time series B is some measure of the artist's "Q" rating. Clearly, increase in A will drive an increase in B, but at the same time, an increase in the Q rating could drive sales higher (or vice-versa). Is there a statistical approach to tackle this problem, or is this more the realm of control systems?


title: correlation between adjacent pixels in an image

text: I would like to know what is the formula for computing the correlation between adjacent values in a 1-D signal. I am reading a paper where there is a formula for converting this to the FWHM of the gaussian distribution and this is given as: $$ (FWHM)^2 = -2 \log(2)/\log(corr) $$ where $corr$ is the correlation between adjacent pixels in the signal. Can I simply do as follows: $$ covar = \frac{1}{N}\Sigma(x_i - \bar{x})(x_{i-1} - \bar{x}) \\ corr = \frac{covar}{\sigma_x * \sigma_{x}} $$ This is the only thing I can think of. Otherwise, I am not sure what it means to compute correlation between adjacent pixels. Can someone verify if this is something that makes sense? Also, why does it make sense? I am not sure why this represents the correlation between two adjacent pixels? Does it even make sense to say correlation of two values???! Thanks for any help you can give me.


title: Coupled hidden Markov model MATLAB error

text: I want to train CHMM using Matlab and I am new to use this software. I am using BNET CHMM which is this:               function bnet = mk_chmm(N, Q, Y, discrete_obs, coupled, CPD)     % MK_CHMM Make a coupled Hidden Markov Model     %     % There are N hidden nodes, each connected to itself and its two nearest neighbors in the next     % slice (apart from the edges, where there is 1 nearest neighbor).     %     % Example: If N = 3, the hidden backbone is as follows, where all arrows point to the right     %     % X1--X2     %   \/      %   /\     % X2--X2     %   \/      %   /\     % X3--X3     %     % Each hidden node has a "private" observed child (not shown).     %     % BNET = MK_CHMM(N, Q, Y)     % Each hidden node is discrete and has Q values.     % Each observed node is a Gaussian vector of length Y.     %     % BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS)     % If discrete_obs = 1, the observations are discrete (values in {1, .., Y}).     %     % BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS, COUPLED)     % If coupled = 0, the chains are not coupled, i.e., we make N parallel HMMs.     %     % BNET = MK_CHMM(N, Q, Y, DISCRETE_OBS, COUPLED, CPDs)     % means use the specified CPD structures instead of creating random params.     %  CPD{i}.CPT, i=1:N specifies the prior     %  CPD{i}.CPT, i=2N+1:3N specifies the transition model     %  CPD{i}.mean, CPD{i}.cov, i=N+1:2N specifies the observation model if Gaussian     %  CPD{i}.CPT, i=N+1:2N if discrete               if nargin < 2, Q = 2; end     if nargin < 3, Y = 1; end     if nargin < 4, discrete_obs = 0; end     if nargin < 5, coupled = 1; end     if nargin < 6, rnd = 1; else rnd = 0; end          ss = N*2;     hnodes = 1:N;     onodes = (1:N)+N;          intra = zeros(ss);     for i=1:N       intra(hnodes(i), onodes(i))=1;     end          inter = zeros(ss);     if coupled       for i=1:N         inter(i, max(i-1,1):min(i+1,N))=1;       end     else       inter(1:N, 1:N) = eye(N);     end            ns = [Q*ones(1,N) Y*ones(1,N)];           eclass1 = [hnodes onodes];     eclass2 = [hnodes+ss onodes];     if discrete_obs       dnodes = 1:ss;     else       dnodes = hnodes;     end     bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'eclass1', eclass1, 'eclass2', eclass2, ...               'observed', onodes);          if rnd       for i=hnodes(:)'         bnet.CPD{i} = tabular_CPD(bnet, i);       end       for i=onodes(:)'         if discrete_obs           bnet.CPD{i} = tabular_CPD(bnet, i);         else           bnet.CPD{i} = gaussian_CPD(bnet, i);         end       end       for i=hnodes(:)'+ss         bnet.CPD{i} = tabular_CPD(bnet, i);       end     else       for i=hnodes(:)'         bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);       end       for i=onodes(:)'         if discrete_obs           bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);         else           bnet.CPD{i} = gaussian_CPD(bnet, i, CPD{i}.mean, CPD{i}.cov);         end       end       for i=hnodes(:)'+ss         bnet.CPD{i} = tabular_CPD(bnet, i, CPD{i}.CPT);       end     end      I put the parameters in mk_chmm(4, 4, 4, 1, 1, 0) in the editor window, but it gives me an error:               Undefined function 'mk_dbn' for input     arguments of type 'double'.          Error in mk_chmm (line 70)     bnet = mk_dbn(intra, inter, ns,     'discrete', dnodes, 'eclass1', eclass1,     'eclass2', eclass2, ..      Any ideas how can I solve it please?? Sorry I am a beginner.


title: Why do different estimators for stock volatility exist? (Realized Variance, RAV, etc)

text: I am very confused about why different volatility estimators (RV, RAV, BPV, etc) exist. If the goal is to find the best estimator for stock volatility, and volatility is latent, how do I know which estimator performs the best when the mathematical definition for these estimators are all different?               Realized variance = sum[(return at t)^2] in a day     Realized absolute value = sum|return at t| in a day     Bipower variation = sum |return at t-1|*|return at t| in a day      If I use a linear model and substitute each of these estimators as the regressor, doesn't this mean I can't compare the results directly, but I can only compare how well the linear model is able to predict itself based on MSE, R^2, and so on?


title: Big Data vs multiple hypothesis testing?

text: Nate Silver in his excellent "The Noise and the Signal" warned that we are much in awe of Big Data. But, that Big Data predictions in many fields have been disastrous (financial markets and economics just to name a few fields). With more data, you get more spurious correlations, more false positives, and erroneous answers. In stating so, he also liens on the excellent work of Loannidis who indicated that over 2/3ds of scientific findings are wrong as they can't be replicated (based on extensive reviews of working papers). In other words, watch out for the many traps of multiple hypothesis testing, especially when you have not even phrased the hypothesis to begin with. "Correlation does not entail causation" still prevails. Now in a new book (called Big Data) written by Viktor Mayer-Schonberger and Kenneth Cukier, Big Data looks far more promising. Given the size of the sample that often equates to the entire population, you can detect granular relationships between subsets of the data you could never before. And, within this Big Data era correlation seems far more important than causation. Figuring out what variables are predictive gets you far better and rich results than figuring out which ones are truly causal (that often turns into an elusive chase). The author mentions several new tools that are aimed at extracting and analyzing Big Data sets including neural networks, artificial intelligence, machine learning, sensitivity analysis among others. Being unfamiliar with any of those (and very familiar with traditional statistics and hypothesis testing in particular), I can't judge if the author statement is accurate (he is not a quant). Do those techniques truly avoid the traps of spurious correlations, multiple hypothesis testing, model overfit and false positive results? Can you reconcile both views: Nate Silver vs Viktor Mayer?


title: Interpreting QQplot - Is there any rule of thumb to decide for non-normality?

text: I have read enough threads on QQplots here to understand that a QQplot can be more informative than other normality tests. However, I am inexperienced with interpreting QQplots. I googled a lot; I found a lot of graphs of non-normal QQplots, but no clear rules on how to interpret them, other than what it seems to be comparison with know distributions plus "gut feeling". I would like to know if you have (or you know of) any rule of thumb to help you decide for non-normality. This question came up when I saw these two graphs: ![graph 2](http://i.imgur.com/UcPyaad.jpg) ![graph 1](http://i.imgur.com/XlPcbjB.jpg) I understand that the decision of non-normality depends on the data and what I want to do with them; however, my question is: generally, when do the observed departures from the straight line constitute enough evidence to make unreasonable the approximation of normality? For what it's worth, the Shapiro-Wilk test failed to reject the hypothesis of non-normality in both cases.


title: Is it valid to select a model based upon AUC?

text: I have plot ROC for several models. These models were used to classify my samples into 2 classes. Using these commands, I can obtain sensitivity vs. specificity plots for each model:               perf <- performance(pred, "sens", "spec")     plot(perf)      Should I rely on the area under the curve (AUC) for each model to conclude which model is better? Other than AUC, should we consider other results so as to conclude which model is better? If yes, how to get AUC with R? Am I right in assuming that "the smaller it is the better is the classification power of the model?"


title: How to do prediction from a linear regression?

text: I am not very good in statistics (ok, I'm really bad), I guess this is a very simple question but I dont understand much of the literature. I have a dataset that is arranged in 2 columns (`var1` is my response variable, `var2` my predictor). Based on Pearson correlation coefficient, I found that my data is correlated (with a significant p-value < 0.01). I also fitted a regression model in R as `lm(var1 ~ var2)`, and made a nice graph. If my two variables are correlated, does this mean that I can use the results from my linear regression to make a prediction, or do I have to rely on another technique? What I am trying to predict is the number of failures by time X in the same software. The regression equation I obtained reads $0.002282\cdot x +0.545751$. The data sample is not very big, it is the cumulative number of failures experienced in a software program through time.               > summary(FAILURES,TOTALTIME)       Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        1.00    4.25    7.50    7.50   10.75   14.00     >cor.test(TOTALTIME,FAILURES)     Pearson's product-moment correlation     data:  TOTALTIME and FAILURES      t = 54.1572, df = 12, p-value = 8.882e-16     alternative hypothesis: true correlation is not equal to 0      95 percent confidence interval:     0.9933656 0.9993741      sample estimates:     cor      0.9979606       Here is a graph of my data: ![scatterplot with fit](http://i.stack.imgur.com/gNZze.jpg)


title: How do you generate ROC curves for leave-one-out cross validation?

text: When performing 5-fold cross-validation (for example), it is typical to compute a separate ROC curve for each of the 5 folds and often times a mean ROC curve with std. dev. shown as curve thickness. However, for LOO cross-validation, where there is only a single test datapoint in each fold, it doesn't seem sensical to compute a ROC "curve" for this single datapoint. I have been taking all of my test data points (along with their separately computed p-values) and pooling them into one large set to compute a single ROC curve, but is this the statistically kosher thing to do? What is the right way to apply ROC analysis when the number of data points in each fold is one (as in the case of LOO cross validation)?


title: Assessing and testing inter-rater agreement with kappa statistic on a set of binary and Likert items?

text: I am trying to calculate inter-rater reliability scores for 10 survey questions (on which there were 2 raters)--seven questions are binary (yes/no), and 3 are Likert-scale questions.   1. Should inter-rater reliability be tested on EACH of the 10 questions, or is there an overall inter-relater reliability test that tests reliability of all questions at one time? If so, what is it?   2. For the binary questions, the agreement level between the two raters is 70-90% on nearly every question, however the Kappa score is often very poor (0.2- 0.4). Can this be right? (And if so, is there a more appropriate test?)   3. Finally, can you use a Kappa-based test on Likert scale questions? If not, what is the correct test for inter-rater reliability?


title: Genetic algorithm for parameter estimation

text: I am a bit confused about parameter estimation using evolutionary methods and their ability to do such a job. Since I am not that pro in stat I am describing my problem with giving an example. Given different response variable, say Y and in my example molecular weight there are different constant value/weight that can lift up our variables (atoms) to hit that response value. The example is molecular weight of molecules. There is 10 molecule which combined by 100 different atom and we are wondering about weight of each atoms. Say $CCOOONH$ is one of the molecules and we know its molecular weight, $\alpha$ so the equation is like $\alpha$= 2$C_w$+3$O_w$+1$N_w$+1$H_w$. My question is manly focus on the relation that GA or other evolutionary methods are able to derive (can they catch a nonlinear relation?!), how much it depends on the fitness function (if I work with a nonlinear fitness function) and finally, could you provide some survey or packages which are relevant? Could you elaborate more about choosing a proper fitness function in your answer, for instance in my question ... could it be a simple correlation function of a new linear combination of variables and the response variable?!   How about ordering variable, meaning some variable could be look like $\alpha$ = 2$C_w$^$\gamma_1$+3$O_w$^$\gamma_2$+1$N_w$^$\gamma_3$+1$H_w$^$\gamma_4$ which apparently has a heavier computation.


title: When do (and don't) confidence intervals and credible intervals coincide?

text: Yes, I know there are many questions on comparing these two types of intervals, but none of them appear to answer this exact question. Here is a blog post demonstrating one case where the two intervals are the same and one where they are different: https://jakevdp.github.io/blog/2014/06/12/frequentism-and- bayesianism-3-confidence-credibility/ But are there any general rules to this? My personal experience suggests that misinterpreting confidence intervals as credible intervals is widespread (I would guess that near 100% are misinterpreted in the medical literature). However, there seems to be no problem with this in practice as long as the confidence interval is a close approximation to the credible interval.


title: Multi-Level Framework for an Evalution Design

text: This is my first time posting here, so please forgive me if I am not following some form of etiquette. My question is in regard to the following evaluation design: We are examining the impact of policy X. It was implemented at the same time for all units j. We have the ability to collect measurements of A for all participants (i) in a given unit j in both the pre- and post-implementation periods. We do not have a non-treated post-implementation group for comparison, so the following was suggested: We can match post-participants to the pre-participants on relevant person level factors. Second, to account for historical changes in j-level factors, we can adjust for changes in said factors at level 2 of a hierarchical model. Overall, our goal would be to estimate an adjusted avearage change in A at the j-level or even estimate what the predicted change would have been when we assume no change in j-level factors. What are the downsides of this approach?.


title: Can I convert a covariance matrix into uncertainties for variables?

text: I have a GPS unit that outputs a noise measurement via covariance matrix $\Sigma$: $\Sigma = \left[\begin{matrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz} \end{matrix}\right] $ (there's also $t$ involved but let's ignore that for a second.) Suppose I want to tell someone else that the accuracy in each direction ($x,y,z$) is some number. $\mu_x, \mu_y, \mu_z$. That is to say, my GPS may give me a reading of $x=\bar{x}\pm\mu_x$, etc. My understanding is that $\mu$ in this case implies that all measurands are independent of each other (i.e., the covariance matrix is diagonal). Moreover, finding the vector error is as simple as as adding errors in quadrature (square root of sum of squares). What happens if my covariance matrix is not diagonal? Is there a simple number $\mu_x^*$ that encompasses the effects of the $y$ and $z$ directions? How can I find that given a covariance matrix?


title: Lognormal Distribution as Maximum Entropy Probability Distribution

text: According to the wikipedia article on the lognormal distribution, the lognormal distribution is "the maximum entropy probability distribution for a random variate $X$ for which the mean and variance of $\log(X)$ is fixed". Is there a not too complicated account of what this means and how this is derived?


title: how to deal with a very predictable variable - but very small sample - should I keep them ?

text: I am doing modeling for say response y ~ a+b+c+e+f.. x+y+z ( say I got model for decision tree ) if my dataset is big, say 1M records : I have a-z independent variables , but having x,y,z are very predictable ( x,y,z are binary variable - 0,1) : say correlation of almost 1. , but very small say ~ 100 customer that would have this variable x,y,z did not show up in the tree, but should I keep them if I know that anyone with x,y,z would have almost 100% response ? Which approach should I use in this case ? Thank you.


title: Algebraic definition of a residual from a regression

text: Can anyone give some advice on how to start proving this algebraically? Define the residual from a regression (one independent variable) algebraically and show that:   1. the mean of the residuals is zero   2. the correlation of the residuals and the independent variable is zero


title: Terminology question concerning regression

text: hopefully you can help me with the meaning of the following, I don't really understand the terminology: "regression of a vector of ones on the matrix $W$", where $W$ is something like $(W_t)' = (w_{1t},w_{2t},w_{3t}, w_{4t})$. I don't understand, which regression I actually have to compute. If it is of help for you, I'm trying to understand the Godfrey/Wickens test according to this article: Ghosh, S., Gilbert, C. L., & Hallett, A. J. H. (1983). _Empirical Economics_ , 8, 63–69. DOI: 10.1007/BF01973190 Greeting, Julius


title: Changes in the regression coefficient

text: "Under what condition (or conditions if you think it necessary) would one observe no change in the regression coefficient (e.g., b-hat Y on X1) for some variable when another variable is added to the regression equation?" I think the answer is when the exogenous variables are perfectly uncorrelated - is that correct?


title: Too many p-value less than 0.0001 is alarming?

text: Can anyone explain why having too many p-value less than 0.001 is alarming? Like what happens to my model right now: ![enter image description here](http://i.stack.imgur.com/rkeDb.jpg) ![enter image description here](http://i.stack.imgur.com/Y1IYY.jpg) All independent variables are categorical except for Calendar_Year. There are 682211 observations in my data set.


title: DCC GARCH in Stata

text: I want to estimate the dynamic conditional correlation GARCH in stata 12. I just ran the following line:               mgarch dcc (x y =, noconstant) , arch(1) garch(1)      The x and the y variables are log first differences from two price series. My questions are as follows:   1. Does anybody know how I obtain the parameter estimates for the elements of Ht?   2. Is it correct that I directly put the log first differences in this model or should I first obtain the residuals from a AR model?


title: Probability of a card being higher than a random card

text: Say you have a game that asks you whether you think the next card will be higher. In this game ace is 1. If the next card is equal to the first one then another random card is taken until the next card is either higher or lower. Suits are disregarded in this game. One deck is used and cards are recycled. What is the probability of the next card being higher than the previous one. I have thought up this equation but not sure if it is correct: $$probability = \frac{13 - Rank}{12} * 100% $$ Where Rank is the rank of the current card (from 1 to 13) Therefore if Rank = Ace, probability = 100% Is this equation correct? Thanks!


title: Example of "real life" use of Bayesian inference on $\mu$ from a normal distribution?

text: A classic example for students, when teaching Bayesian statistics, is to make inference on the mean parameter $\mu$ of a normal distribution, when it has a prior normal distribution. I would like to find a real-life example of why such a construction would be not only elegant, but also useful. Can anyone direct me to such an example?


title: GMM weight matrix $W_n$ -- what does the index $n$ signify?

text: The question is as in the title. The GMM estimator (in my Econometrics notes) is as follows: $\hat \delta_{GMM}(W_n) = (S'_{xz} W_nS'_{xz})^{-1}(S'_{xz}W_ns_{xy})$


title: IVEware warnings

text: When using SAS callable IVEware for multiple imputation it will occasionally throw out a warning for too many iterations. Can someone give me an idea of how many warnings are acceptable (if any) for a particular variable?