title: Structuring a linear mixed model in R with nesting

text: My ecological question is: "What are the trends in percent coral cover by island and depth across the state of Hawaii from 1999 to 2012?" I am trying to analyze this hierarchical data set using R with 10 transects at each depth, 2 depths per site, and site nested in island. Data structure:               Fixed effects:      Island: Hawaii, Maui, Molokai, Kahoolawe, Oahu, Kauai.      DepthCat: S = Shallow, D = Deep.      WYear: 0-13. It was suggested that I use this factor as a covariate for years.          Random effects:      Site: 34 sites across the 6 islands with 2 depths per site.      Transect: 10 permanent transects at each depth.      Year: 1999 – 2012 (14 years)          Dependent variable: PercentCover      Currently, I am using the `lmer` function in the `lmerTest` package and this is the model that I've constructed.               fit1 <- lmer(PercentCover ~ WYear*Island*DepthCat +                  (1+WYear|Island/Site/DepthCat/Transect) + (1|Year), data=Benthic)      Unfortunately, the data are spotty (i.e., missing data in multiple years for a number of sites) so the model returns `[1] "Asymptotic covariance matrix A is not positive!"`, even using arcsin transformed data. I can still run the summary statistics to get results, but I don't feel comfortable with the error message. Perhaps I have not structured the model correctly in terms of organizing the nested factors, but the number of observations for each of the levels in the summary stats seems correct. I tried different and simpler iterations of the model such as:                fit1 <- lmer(PercentCover ~ WYear + Island + DepthCat + (1+WYear|Transect/Site) +                    (1|Year), data=Benthic)      which works, but doesn't give me the interaction information and returns a larger AIC suggesting that the model does not fit the data as well. To deal with all of the missing data, I tried another approach by using the regression slope of percent cover over time as the dependent variable for each site X depth combination.               Data structure:          Fixed effects:      Island: Hawaii, Maui, Molokai, Kahoolawe, Oahu, Kauai.      DepthCat: S = Shallow, D = Deep.          Random effects:      Site: 34 sites across the 6 islands with 2 depths per site.      Transect: 10 permanent transects at each depth.          Dependent variable: Trend      I used the following model, but the summary results did not make much sense, even after transforming the data.                fit1<-lmer(Trend ~ Island*DepthCat + (1| Island/Site/DepthCat/Transect), data=Benthic)      Any suggestions on improving my analytical approach would be appreciated.


title: Two-sample location test when data is heavily skewed towards zero?

text: I want to test if the mean between two independent samples is different. Both samples are large, about 2 million observations each, however almost all of the observations are zero. In particular in each sample, only around 10,000 observations are non-zero. This data captures the amount of money spent over a given time period amongst two groups of customers, a value of zero means that the customer did not make any purchases. Running a Student's t-test or Mann-Whitney U test does not show any statistically significant difference between the means. I was considering dropping all of the non-zero observations, and re-running the tests, but I'm not sure how to interpret the results. Any insight would be appreciated.


title: Behavioural economics experiment - non parametric test

text: I wonder if somebody could help with my dissertation project. I am stuck with the analysis and tried to consult a few people but unfortunately they either have limited experience, knowledge or interest in the subject of statistics/behavioural economics. Basically, I got a mixed repeated measures design with 2 dependent variable (DV) and 2 independent variable (IV) (1 at 2 and 1 at 4 levels). There were 24 participants split between two groups. Each participant performed a task which was measured using two DV between groups. Also, in each group there were four levels of additional IV. Given that I must use non-parametric tests (non-normal distribution) and i have 2 DV (and so MANOVA) I do not think I could use the SPSS for this. I am not sure which tests to use and what software (R?).


title: Fitting a back-to-back distribution

text: This may be a very simple question, but what does it mean "to fit a back-to- back Weibull distribution" to a residual series (or I suppose, to any arbitrary data)? I think it means to fit a series like this: $x_{t} = \sum_{k}\alpha_{k}\frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{\left(\frac{x}{\lambda}\right)^{k}}$ but I don't really know. **Context:** From a problem in Cowpertwait's book "Introductory Time Series with R": > A hydrologist wishes to simulate monthly inflows to the Font Reservoir over > the next 10-year period. Use the data in Font.dat (§2.3.3) to answer the > following: > > Plot a histogram of the residual errors of the fitted AR model, and comment > on the plot. Fit a back-to-back Weibull distribution to the residual errors.


title: How to filter these extremely bad data points?

text: I'm looking at data for my company, and basically we have some periods over the last year where the data was not uploaded correctly. In this figure, "mu" is the value of interest, and duration is a time in days. The huge spike near 180 days is an artifact, as are the smaller spikes at (140, 200, 275, etc). Because I don't want to smooth out the initial part of the curve near 0 days, I'm having trouble coming up with a good way of smoothing out these artifacts, and have come here looking for suggestions. Thanks, in advance! ![enter image description here](http://i.stack.imgur.com/8z1Kj.png)


title: Bayesian Classification evaluation

text: I am trying to implement Bayesian Classification on the data set as follows: "Problem: classify whether a given person is a male or a female based on the measured features. The features include height, weight, and foot size." Example training set below (borrowed from http://en.wikipedia.org/wiki/Naive_Bayes_classifier) ![enter image description here](http://i.stack.imgur.com/3fprY.png) **Testing** Below is some samples to be classified as a male or female.               sex       height(feet)  weight(lbs)      foot size(inches)     sample         6            130            8     sample         4            140            5      **Results:**               Actual:'female' and the output is (0.99, 'female')     Actual:'female' and the output is (0.89, 'male')      Here, i want to evaluate the classification accuracy. I am confused between two options for evaluation: First option, it classified 1 of 2 test data correctly, so the accuracy is 50%. Second option, the probabilities of being female are (0.99 + 0.11) /2 = 0.55, so the accuracy is 55%. Which option makes sense to evaluate the classification accuracy of this method? What are the meaning of each options?


title: Automatic outlier detection in R

text: Our model processes millions of multivariate observations; manual outlier detection is impractical. I am looking for a method of automatic outlier detection. I have been trying to use R package `mvoutliers`, especially function `pcout`, and get the error > More than 50% equal values in one or more variables! The problem is that our data is quite sparse and many variables include _More than 50% equal values_.   1. Is there any way around (some data preprocessing) that would still allow me to use `pcout`?   2. Is there another recommended R package/function/method for automatic outlier detection?


title: "Constant variance" violation

text: If we apply linear regression on a data which has a BINARY(0,1) dependant variable, the very important assumption of "constant variance" of the dependant variable across independant variables is violated. Can anyone explain how ?


title: How do I find values not given in (/interpolate in) statistical tables?

text: Often people use programs to obtain p-values, but sometimes - for whatever reason - it may be necessary to obtain a critical value from a set of tables. Given a statistical table with a limited number of significance levels, and a limited number of degrees of freedom, how do I obtain approximate critical values at other significance levels or degrees of freedom (such as with $t$, chi-square, or $F$ tables)?


title: Measurement of similarity for hierarchical clustering trees

text: I am performing a number of hierarchical clusters on a dataframe of patient records (e.g. similar to http://www.biomedcentral.com/1471-2105/5/126/figure/F1?highres=y) I am experimenting with different distance measures and different parameter weights, to understand the impact on the final clusters / structure of the tree (dendrogram). My question whether there is a standard calculation / measure to calculate the difference between different hierarchical trees, and how to implement this in R (e.g. to quantify that some trees are nearly identical, and that some are drastically different).


title: How do I go from Labor Statistics results to generalities on the population?

text: I'm currently looking at survey results from the US bureau of labor statistics concerning hourly wages. With the raw survey data, I get different respondents who state how many hours they work in a week, along with the average pay, and some demographic data like age and sex. How can I go from this information to making an educated guess about the amount of hours have been worked by the entire US population at a certain wage level, for example? Is there some way to "correlate" the values considering more than just a single parameter like age? EDIT: I'm going to try and ask a more specific question then. Please forgive my complete lack of statistics knowledge: Let's say I'm doing a survey about how much people earn in an hour. Based off of my survey, I know that , in my survey sample, men make about \$10 an hour, and women make \$8 an hour. I also know that in the same sample, people under 40 make \$7 an hour on average, and people over 40 make \$11 an hour. Assuming that the population is 50/50 male/female, I can make an educated guess that the population , on average, makes \$9 an hour. Likewise, if %30 of the population is over 40, then I can guess that the population on average makes \$8.20 an hours But how can I "combine" both of these information points? How can I consider both the sex and the age of the people to fit it to the population? Am I just supposed to do a weighted average or something?


title: Multiple dependent variables with categorical predictors: estimation problem

text: I have two continuous variables which are my dependent variables (well they are discrete actually) and two predictors which are categorical (they are experimental treatments). `O + M ~ P + Q` I have already broken the model into pieces and I have tested through two-way factorial ANOVAs the effect of the predictors P + Q on each of the dependent variables, O and M, separately: `O ~ P + Q + P*Q` and `M ~ P + Q + P*Q` The dependent variables are two attitudinal items on different issues (they are both agreement to a certain statement, statement O and statement M, which are very different with each other). P, the first predictor, is a classification of belief system derived from a self-reported measure (namely the person self-set him/herself as conservative or liberal) it is a bit endogenous with O and P but not always the case (here I treat it as it is not). The variable Q is the most important one because it's the experimental treatment (split-sample); it gives us information on how the data where collected (context). By deciding to interact P and Q I presuppose that in specific type of Q each level of P will be affected differently in terms of the outcome variables (O and M). This was the case for attitude O, but not for the attitude M. Now I would like the combine the two samples which were divided by Q (the split-sample experimental treatment) and see whether P and Q interact when on both variables at the same time. And if yes, at how my response variables are affected. I want to examine the effect on both response variables O + M which I already know that they are negatively correlated when predicted by P, namely, the measurements in O are likely to be the opposite than in M in all levels of P (I mean the means). As I said before I have alredy tested it in two-way factorial ANOVAs and it was significant. Moreover, something rare occurs (I have an explanation from my field) the Q predictor/experimental-treament is significant for the M dependent variable but not the interaction term P*Q Due to the needs of my project I need to run the model mentioned above `O + M ~ P + Q` but also with the interaction term `O + M ~ P + Q + P*Q`. My first question is what method should I use? I have already ran this in R, but it truncates my predictor in the `summary.lm()` output, namely, it gives me estimates only for 2 from the 7 categories of P... From what I've read I need to apply a MANCOVA, however, is it still called a MANCOVA if there are two categorical predictors and their interaction term included? Also, what is the type of analysis (term) for getting estimates on my dependent variables (O + M) for each level of P and Q and their interaction term? What is this type of analysis? Can you suggest any R package? Thank you for your help


title: Conditional Probability with Normal Distributions?

text: Let's say that I have $3$ independent random normal variables, $A$, $B$ and $C$. They all have a standard deviation of $17.526$, while $A$ has a mean of $143$, $B$ of $139$, and $C$ of $129$. My question is, how would I go about calculating $P(A>C\mid A>B)$? I know how to calculate $P(A>C)$ and $P(A>B)$, but it's the conditional part I'm struggling with. I'm not sure how I would manage to apply Bayes' theorem to this, so I'm stuck. If it helps speed-wise, I have $P(A>C)$ to be $0.713912$, and $P(A>B)$ to be $0.564105$ with the difference distribution in each case having a stdev of $24.785$. Thanks.


title: Weibull regression with known intercept in R

text: I'd like to run a Weibull regression with the pre-defined scale and shape parameters of the Weibull distribution. I'm using `survreg()` from `library(survival)` in R and as I understand in the `survreg` terminology:               scale = 1/(rweibull shape)     intercept = log(rweibull scale)      Specifying the `shape` parameter (sic!) by running `survreg(Surv(y)~x, scale=1/rweibull_shape)` seems to work and, as expected, affects mainly the significance of the regression coefficient for the explanatory variable and not the coefficient itself. As for the scale parameter I assumed I needed to use `offset`:               survreg(Surv(y)~x+offset(rep(log(rweibull_scale),length(x))), scale=1/rweibull_shape)      However, specifying the scale this way doesn't seem to affect the model's log likelihood and the significance of regression coefficients, which is surprising given that Weibull variance depends on both shape and scale. What am I doing wrong? Full example:               ## offset –0.1     summary(survreg(y)~x+offset(rep(-0.1,length(x))), dist="weibull", scale=0.143))          Call:     survreg(formula = Surv(y) ~ x +      offset(rep(-0.1, 4)), dist = "weibull", scale = 0.143)                        Value Std. Error    z        p     (Intercept)        0.690     0.1599 4.31 1.61e-05               x        0.406     0.0715 5.68 1.31e-08          Scale fixed at 0.143           Weibull distribution     Loglik(model)= -6.5   Loglik(intercept only)= -15.1     Chisq= 17.22 on 1 degrees of freedom, p= 3.3e-05      Number of Newton-Raphson Iterations: 7      n= 4                ## offset –10     summary(survreg(y)~x+offset(rep(-10,length(x))), dist="weibull", scale=0.143))     Call:     survreg(formula = Surv(y) ~ x +      offset(rep(-10, 4)), dist = "weibull", scale = 0.143)                             Value Std. Error     z        p     (Intercept)        10.590    0.1599 66.23 0.00e+00               x        0.406     0.0715 5.68 1.31e-08          Scale fixed at 0.143           Weibull distribution     Loglik(model)= -6.5   Loglik(intercept only)= -15.1     Chisq= 17.22 on 1 degrees of freedom, p= 3.3e-05      Number of Newton-Raphson Iterations: 7      n= 4


title: PC analysis and regression-percentage of explained variance

text: my question is about the percentage of explained variance in principle component. different codes in R software show different values for percentage of explained variance(PEV) with the obtained PEV by hand.(the PEV for each principle component by hand is:73 % 25% the formula to calculate PEV is: eigenvalue/sum(eigenvalues) for example using the following codes(prcomp, pcr) on the convariance matrix S v1 v2 v3 [1,] 1 -2 0 [2,] -2 5 0 [3,] 0 0 2 > summary(prcomp(S))$importance[3,] PC1 PC2 PC3 0.92438 1.00000 1.00000 > > # using pcr code > > y<-c(2,3,4) pcr(y~S) explvar(mpcr) Comp 1 Comp 2 92.43819 7.56181 Thanks


title: modelling data that contains only ones and zeros

text: I'm new to modelling this type of data, and I got punished for asking this on stackoverflow... I have a dataset where the predictive variables contain only ones and zeros, and the response variable also contains only ones and zeros. I would like to build a model in R to describe their relationship. I know that a logistic model is suitable for a binary response variable, but what is an appropriate model for binary predictor variables?


title: Binary outcome with categorical predictors - experiment control

text: I have a problem which I try to formalize/decompose/map to well-known cases, so this is more of a call for references/names/terms. Description of the problem: The system receives a stream of plant seeds, each seed having a few categorical properties. For each seed, the system chooses a cultivation protocol to use (thus adding a few more categorical properties). As a result of cultivation, a seed may eventually bring a single fruit (but more likely will never bring any). The ultimate goal of the system is to maximize the number of fruits. End of description. Ideally, I would like a pointer to literature concerned with a similar problem (not necessarily about plants, may be formulated for manufacture similarly). I recognize, though, that the problem is a complex one, and it is more likely to find prior art for separate pieces of it. There are also some simplifications I am willing to allow. E.g., instead of any seed being in either "fruit is obtained" or "fruit is not yet obtained but potentially may be" states, it can be stated as being in either "fruit is obtained" or "fruit is not obtained and never will be", thus reducing cultivation of every seed to a Bernoulli trial. Let's call the probability of getting a (single) fruit "fecundity" (though this is different from the usual meaning in biology). One of the possible decompositions I see is reducing the goal of maximizing total number of fruits to density estimation of fecundity (indexed by all possible combinations of properties of both seed and its cultivation protocol). This does not answer how to combine exploration with exploitation, but I hope to use the multi-armed bandit analysis for that. An extremely simplified example would be 2 properties for seeds: color (red, black, white), pattern (striped, solid, spotted) and 1 property for cultivation: soil (sandy, mud, hydroponic). Upon seeing a red striped seed the system could decide to plant it into sandy soil, for example (or maybe to plant it into sandy soil 60% of cases, and 40% into mud). In real problem the number of properties is ~ 6-7, and some of them have hundreds/thousands of values (e.g., source farm). This is already bad, but in addition, fecundity is very low (~1 fruit per 10 thousand seeds), so even with millions of experiments, the coincidence matrix of property values is going to be quite sparse - thus I am not optimistic about using linear logistic regression. I think the only way to arrive to any estimates of fecundity for most of combinations of seed/cultivation properties is by utilizing dependency among fecundity for "similar" combinations - by similar here I mean the ones sharing at least one property value. Empirically, existence of such dependency is confirmed, but its specific nature is unknown - e.g., fecundity of seeds of the same color on different cultivation protocols more closely correlates to each other than to fecundity of random seeds; and the same is true for seeds of the same pattern, but it is not clear how to model this (hierarchical Bayesian generative model?). At the moment I am looking at http://en.wikipedia.org/wiki/Thompson_sampling, with a context corresponding to a seed, and an action corresponding to a cultivation protocol, but the main missing piece is the specific model of distributions. Any ideas or even questions are appreciated!


title: Linear model for data that follow gaussian distribution

text: I have a question about linear regression. We have the linear regression of input data $(X,Y)=((x_1,y_1),(x_2,y_2)...(x_n,y_n))$ is $$F=aX+b$$ a,b are factors of the linear line, $y_i$ is {-1,1}. Given the probability ofx is normalization distribution $$p(X|Y=1 ;-1,a,b) is N(\mu,\sigma)$$ How to find a,b satisfy that $$L=\sum(Y-F)^2$$ is minimum Thank you so much. Could you suggest to me the paper about that issue


title: the combination of two independent continuous time Markov chains

text: Assume one variable $x$ has two states 0 and 1, $x$ changes between 0 and 1 following a continuous time Markov chain. The transition probability is represented as matrix $P$ and the time sojourning on a state follows exponential distribution with arrival ratio $\lambda$. There is the other variable y which is independent with $x$. $y$ also changes between 0 and 1, and have the same $P$ and $\lambda$ with $x$. Now the question is: will the new variable $z=x*y$ follow continuous time Markov chain? if yes, what will be the transition probability and $\lambda$. ($z$ only has two states: 0 and 1.) I guess $z$ still follows continuous time Markov chain, but I have no idea how to prove? (maybe there has been a property stating this, which I don't know?) Could anyone give some ideas? Many thanks!


title: Negative Weak Correlations

text: I've performed correlation tests and received negative weak correlations using ordinal and interval data. I've tired a log transformation, and still the correlations were the same. I'm looking for a positive weak to moderate correlation with the data. So far, I've done a log transformation, Spearman and Pearson Correlation. You assistance is much appreciated. Thank you!


title: Determining optimum sample size and statistical tests

text: I'ma stats newbie and do not know how to go about this. So I will lay it out as best as i can in the hope that some one here can help me out ! I work for a IT Infrastructure (Datacenter services) company XYZ. This company has sets up wintel, Unix boxes, Databases and Networks for it's clients as part of its services offering. I am trying to determine whether we can do some statistical test using sampling to test the quality of the the work performed by our engineers(Hardware setup plus software installation+any software programs written as part of all this work). So let's say we provide 5 different services (wintel, Unix, DB, Networks, Middleware). All these have two kinds of projects. Standardized projects(like a shrink-wrapped product) and Non-standard/custom projects. Please let me know if you need any more details. If I can use statistical tools for testing, what kind of tests would they be ? what kind of sampling mechanism can I use to determine an optimum sample size ? The idea is to use these tests, get results, look at root causes and then maybe run some six sigma type improvement programs for the deficiencies. Thanks in advance to your replies ### Edit (in response to comments) Broadly we would like to measure "Quality of Performed Work" We are thinking of two types of quality indicators. 1) One would be done at the end of the project - customer survey. 10 questions; responses(5-1 points for Excellent, very good et al) 2) One quality indicator to check during progress of project. for a 3 month project we could have a safety checklist of 4 items, one safety check per every 3 weeks. The idea is to compare is "in-process" project safety check scores with the customer survey scores to get an integrated view of the course of the project. This would help identify if and where a problem has occurred, customer perception, PM communication, engineer execution, process stability, insufficient safety checks etc. Since this entire process is in conceptual stage, there is only limited information about it with me. I cannot characterize it's shape right now. I am assuming that Let me see if I can get some information on that. ### Edit (migrated from a reply that functioned as a comment) Michelle, The Org. i work for is a Matrix organization with the a director level person serving as the Project Sponsor. We are also heavily ITIL based organization where everything is a Service. The Engineers have their Direct reporting Managers and also have dotted line reporting to the Project managers. In this situation where all org.s have their priorities, it will be very hard to tie at Engineer level , their performance goals to some Quality measures. But the VP has defined some Goals for the org, part of which is that everything should be measurable and Imp. measures should be on a Dashboard. So compliance to this will be measured for each group under the org. **I would like some more advice on the kind of qualitative measures that we can use in this approach (as stated in my last two posts). This is to measure Project Delivery by the org.** Thanks to all who contribute. Rolando2- Thank you for your contribution !


title: Hidden Markov models with Baum-Welch algorithm using python

text: I'm looking for some python implementation (in pure python or wrapping existing stuffs) of HMM and Baum-Welch. Some ideas? I've just searched in google and I've found really poor material with respect to other machine learning techniques. Why?


title: Non-exclusive categorical model

text: I have a project to model users with characteristic tags (e.g. runner, cyclist, swimmer, vegan, pianist) in order to correlate user behaviour to these labels. Obviously a user can have multiple characteristics (non- exclusive) and so the distance (or similarity) between users is determined by the amount of tag overlap.   1. Is there a name for this kind of model?   2. Are there any particular patterns/algorithms that are useful for measuring correlations in data represented by this kind of model? I don't need a complete list; just to be pointed in the right direction.


title: Creating joint conditional probability distribution

text: I have a signal gathering system. Every 10 seconds it goes to locations A, B and C and sees whether there is activity or not. It then populates a table which looks like following:               A    B    C     ------------     1    0    1     1    1    0     1    0    1     1    0    1     ............      There are about 10,000 such rows of data in table. Now I want to understand $P(B,C|A)$. My understanding is that I would have 8 cases here like the following: $P(B=1,C=1|A=1)$, $P(B=1,C=1|A=0)$, $P(B=1,C=0|A=1)$, etc. My questions are :   1. Should I model this as binomial distribution or Dirichlet distribution. Why choose one over the other for this case?   2. I know that the way logic has been coded, B and C are impacted by activity at A. But what if C is also impacted by activity at B. Does considering binomial or Dirichlet distribution factor any interaction between B and C in calculating $P(B=1,C=0|A=1)$ for example? I am pretty much of novice in this field so please explain your answer. This question came out of the discussion on this previous question: Creating conditional probability distribution


title: Can mutual information gain value be greater than 1

text: I have a very basic doubt. Sorry if this irritates few. I know that Mutual Information value should be greater than 0, but should it be less than 1 ? Is it bounded by any upper value ? Thanks, Amit.


title: When to choose PCA vs. LSA/LSI

text: **Question:** Are there any general guidelines with respect to the input data characteristics, that can be used to decide between applying PCA versus LSA/LSI? **Brief summary of PCA vs. LSA/LSI:** Principle Component Analysis (PCA) and Latent Semantic Analysis (LSA) or Latent Semantic Indexing (LSI) are similar in the sense that all of them rely fundamentally on the application of the Singular Value Decomposition (SVD) to a matrix. LSA and LSI are, as far as I can tell, the same thing. LSA differs from PCA not fundamentally, but in terms of the way the matrix entries are pre- processed prior to applying the SVD. In LSA the preprocessing step typically involves normalizing a count matrix where columns correspond to 'documents' and rows corresponds to some kind of word. Entries can be thought of as some kind of (normalized) word-occurrence- for-document count. In PCA the preprocessing step involves computing the covariance matrix from the original matrix. The original matrix is conceptually more 'general' in nature than in the case of LSA. Where PCA is concerned the columns are usually said to refer to generic sample vectors and the rows are said to refer to individual variables that are being measured. The covariance matrix is by definition square and symetric and in fact it is not necessary to apply the SVD, because the covariance matrix can be decomposed via diagonalization. Notably, the PCA matrix will almost certainly be denser than the LSA/LSI variant - zero entries will only occur where the covariance between to variables is zero, that is where the variables are independent. Finally one more descriptive point that is made fairly frequently to distinguish the two is that > LSA seeks for the best linear subspace in Frobenius norm, while PCA aims for > the best affine linear subspace. In any case, the differences and similarities of these techniques have been hotly debated in various forums throughout the internets, and clearly there are some salient differences, and clearly these two techniques will produce different results. Thus I repeat my question: Are there any general guidelines with respect to the input data characteristics, that can be used to decide between applying PCA versus LSA/LSI? If I have something resembling a term-document matrix will LSA/LSI always be the best choice? Might expect to get better results in some cases by preparing the term/doc matrix for LSA/LSI and then applying PCA to the result, instead of applying the SVD directly?


title: R's randomForest can not handle more than 32 levels. What is workaround?

text: R's randomForest package can not handle factor with more than 32 levels. When it is given more than 32 levels, it emits an error message: > Can not handle categorical predictors with more than 32 categories. But the data I have has several factors. Some of them have 1000+ levels and some of them have 100+. It even has 'state' of united states which is 52. So, here's my question.   1. Why is there such limitation? randomForest refuse to run even for the simple case.              > d <- data.frame(x=factor(1:50), y=1:50)     > randomForest(y ~ x, data=d)       Error in randomForest.default(m, y, ...) :        Can not handle categorical predictors with more than 32 categories.      If it is simply due to memory limitation, how can scikit learn's randomForeestRegressor run with more than 32 levels?   2. What is the best way to handle this problem? Suppose that I have X1, X2, ..., X50 independent variables and Y is dependent variable. And suppose that X1, X2 and X3 has more than 32 levels. What should I do? What I'm thinking of is running clustering algorithm for each of X1, X2 and X3 where distance is defined as difference in Y. I'll run three clusterings as there are three problematic variables. And in each clustering, I wish I can find similar levels. And I'll merge them. How does this sound?


title: Understanding biplot for compositional data

text: 1. Say in the biplot of some compositional data, I discovered that the links $\vec{AB}$ and $\vec{CD}$ intersect at roughly 90 degrees. So I can say that $\log {A\over B}$ and $\log{C\over D}$ are roughly uncorrelated. The questions is... so what? What's the use in discovering this?   2. Now, say in the biplot I discovered vertices E, F, G are collinear. Then my 'textbook' says they form a 1-d biplot. So they have unidimensional variability. The question is... what does this even mean?


title: Variance of forecast's error

text: Is there any situation in which $$\text{var}(d_t-f_t)=\text{var}(d_t)+\text{var}(f_t)-\text{cov}(d_t,f_t)$$ can be used instead of $$\text{var}(dt-ft)=\text{var}(d_t)+\text{var}(f_t)-2\text{cov}(d_t,f_t)$$ where $d_t$ is a times series and $f_t$ is it's forecast. thank you.


title: Why is the anova p-value on these series of numbers not significant

text: I have three different samples (three replicates per sample) and lots of variables (1000+) This variable in particular had the largest VIP number from an Orthogonal signal correction PLS-DA of these three different samples. Usually I log2 transform these, but for this example I did not (it doesn't seem to make a big difference)               A<- c(7.98E8,4.88E8,1.68E8,2.21E7,5.37E7,1.41E7,4.21E8,8.03E7,3.67E6)     condition<-factor(c(0,0,0,1,1,1,2,2,2))     variable<- data.frame(raw=A,condition=condition)      when I look at the means they seem really different (especially 0 vs 1) It’s almost an order of magnitude               tapply(variable$raw,variable$condition,mean)                  0         1         2      484666667  29966667 168323333       and the box plot               boxplot(variable$raw ~variable$condition)      ![box plot](http://i.stack.imgur.com/ZqqOR.png) when I try a oneway test               oneway.test(variable$raw~variable$condition)          data:  variable$raw and variable$condition     F = 2.945, num df = 2.000, denom df = 2.702, p-value = 0.2095      the p-value is not even close to being significant. Why is this? Is it the large variance, small number of samples? Maybe there is something I don’t understand about anova because the means seem really different across these samples. I do check that the data seems normally distributed so a parametric test seems reasonable. Maybe there is a different test I should be using?


title: Hamming distance for strings with different length

text: I am working on a project that involves computing similarity metrics between strings. I was would like to know whether it is possible to use hamming distance on strings with differnt length, and if possible, how to go about it. I step by step explanation would be grately appreciated. Thanks.


title: Selecting a Bayes shrinkage prior

text: I'm looking for a way to integrate prior knowledge about a parameter in a context equivalent to Bayesian hierarchical models. I come from frequentist background and I'm uninitiated in hierarchical models, but I've got a gut feeling that this could solve my problem. Given an $L_1$ shrunken OLS regression problem I assume that the regression coefficients are normally distributed $$ \beta_1 \sim {\mathcal N}(\mu,1), $$ and its mean follows a double exponential (Laplacian) distribution $$ \mu \sim {\rm}(\frac{1}{2c} exp{(-\frac{|\mu|}{c}})) $$ Now, suppose I know from previous studies that a certain effect size has been reported. How can I integrate that knowledge into choosing the correct parameter $c$ for the Laplacian distribution? Any thoughts or ideas? I'm using this in the context of LASSO regression so I don't see how I could easily specify a hyperprior for $c$


title: Alternative to One-way ANOVA unequal variance

text: I would like to compare the means across three groups of equal sizes (equal sample size is small, 21). The means of each group _are_ normally distributed, but their variances are unequal (tested via Levene's). Is a transformation the best route in this situation? Should I consider anything else first?


title: How to calculate threshold level for mutual information scores?

text: Im quite new to the world of statistics, and only just finished work on my first (successful) attempt to calculate mutual information scores. Now lets say that for the column pairs AB, BC, CD and DE I have mutual information scores 0.1111, 0,9999, 1.23 and 1.5. I need to use these values to accertain which of column pairs AB, BC, CD and DE are dependent of each other, and which are not. To do so, I need to colculate a threshold limit. I assume that if the threshold level is X, the all MI scores larger than X means that the two columns are related, while scores < X means that the two columns are not related. Can someone point me towards a suitable threhold calculation process?


title: Robust variance/covariance matrix in Poisson regression

text: Suppose I have survival data with more than one row per subject, because I have splitted the follow-up time of each subject into pieces (maybe because I have one or more time-varying variables or maybe just because I want to fit a Poisson model with a non-constant hazard over time). Do I have to use the robust variance / covariance estimator (that is implemented for example in Stata with the option `vce(cluster clustvar)`) to take into account that I have more than one observation per subject (i.e. they're not independent)? **Edit (15 March 2012)** : Lambert and Royston in their book carry out this analysis: they split each subject's follow-up on one time-scale (let's say attained age)[*] and fit a Poisson regression including attained age as the dependent variable (modelled for example using splines) plus the offset, so that it's possible to model the incidence of some disease according to attained age. They **do not** use the robust variance/covariance estimator, but I haven't found in the text any explanation of why the single rows (or _episodes_ ) can be considered independent. **The question:** Can someone explain to me why the single rows (or _episodes_ ) can be considered as independent? [*]To clarify what's been done, let's take for example subject number `1001`. He/she enters the study at `80.00219` years of age and develop the disease at `85.037236` years (`_d==1`). This is what happens to the record of this subject after splitting it up. (The offset variable is defined as `ln(_t- _t0)`)                 id   _d          _t        _t0       1001    0          81   80.00219       1001    0          82         81       1001    0          83         82       1001    0          84         83       1001    0          85         84       1001    1   85.037236         85


title: Using secondary data require normality test?

text: I am undergraduate student and I have a problem in my research. Recently, my supervisor told me that if using secondary data no need to perform normality test, why? It is got any citation or previous research prove this?


title: Using autocorrelation plots to choose the number of inputs for a neural network predicting time series

text: A neural network applied to time series needs to have the number of input nodes defined. Each input is applied to a time point previous to the current point being predicted. If $D$ is the number of past time points taken, how is this number determined? Can the autocorrelation function be used for the lag correlation plot to give insight as to the best number of previous values to use? If so how?


title: Outliers spotting in time series analysis, should I pre-process data or not?

text: My question builds on a previous post on outlier detection in generic time series, and specifically on the answer provided by the always great Rob H. I work for a small-sized manufacturing company that currently handles the issue, i.e. detecting outliers in sales data time series, employing a (dubious) automated off-the-shelf software procedure. I think this kind of approach is questionable at best and, more often than not, I'm not happy with the results I get. I would therefore like to "double check" the output from our software using some alternative method. Rob's idea seemed reasonable, straightforward and easy to implement, so I decided to give it a try. Question is: what if my time series are not "generic"? Stl decomposition highlights a strong seasonality and a varying trend in my data: ![enter image description here](http://i.stack.imgur.com/fGuz7.jpg) (BTW I used `stl(x,s.window="periodic")` like Rob suggested, but IMHO `stl(x,s.window="periodic",robust=TRUE)` would be a better choice since outlier detection is the issue at hand here. Also I'm not really sure about the `s.window="periodic"` part, I tried experimenting with different values a bit, but I don't know how to interpret results. Maybe someone can point me in the right direction?). Back to my question, mine being sales data, the seasonal pattern is (or I think it is) strongly affected by calendar effects. Also I have reason to believe the big level shift in 2009 is due to the financial crisis and it has nothing to do with trend. What do I do here? Should I let the model handle this, or should I pre-process data? Do I perform working-day adjustment and re-allign (is there such a thing?) before-2009 and after-2009 data, or do I let STL decomposition do the work? I could write another 1000 lines, but I think this should be enough to get the message through. I apologize for the WOT and for my bad english. Also I hope I did not break too many forum rules... I hope someone out there can help!


title: Estimate fatigue life parameters in R

text: I have project where I'm given a set of times between a certain event and then I have to estimate the parameters for the ditribution function with the best fit using R, similar to what EasyFit does. Unfortunately, I'm fairly new to R and statistics in general so I'm having a hard time getting off the ground. I'm attempting to use the fitdistr package but I can't figure out if it supports fatigue life functions. Is there another package I'm missing that might serve better? Any help you can provide would be greatly appreciated. Thanks, Nick


title: Survey-Based Analysis, Small Sample Size, Ordinal Dependent Variable

text: I conducted a survey and only received 50 responses. My dependent variable is ordinal. My independent variables are ordinal, dummy, and continuous... What statistical model should I use to analyze the results?


title: Simple similarity metric

text: I am trying to find a simple similarity metric that will compare vectors of indeterminate length. These vectors will be populated with values between 1 and 5. In this situation a 1 is closer to a 2 then it is to a 5 etc etc. I am new to this type of math. I naively considered using cosine similarity. However, when I looked into this further I realized that this probably wasn't the right metric for this sort of calculation. Thanks in advance and I apologize for the newbie question. Note: I will be programming this in PHP.


title: How to use Stata to pool Cohen's d?

text: I am attempting to calculate Cohen's d and then pool those estimates into a summary effect size. Can anyone help? (Stata or SPSS software owned).


title: On the use of oblique rotation after PCA

text: Several statistical packages, such as SAS, SPSS, and R, allow you to perform some kind of factor rotation following a PCA.   1. Why is a rotation necessary after a PCA?   2. Why would you apply an oblique rotation after a PCA given that the aim of PCA is to produce orthogonal dimensions?


title: Using either AIC or BIC to compute regularization parameters

text: Could someone out here kindly help me to state whether or not I am wrong on how to compute the tuning parameter for penalized neural network training problem... Given the objective function, I performed $N=100$ runs, and for each chosen tuning parameter, I obtained an optimal set of parameters (which gives me the required model). Using each set of model, I constructed the following: $\text{AIC} = \text{loss} + \text{complexity}$ and then proceed to determine such a minimum having obtained the entire values for $\text{AIC}$ from the $N=100$ runs. My question is: Is this similar to minimizing using grid search? I now assume the minimum gives me a candidate value for the optimization parameter.


title: How to efficiently generate Positive Semi Definite correlation matrices?

text: I would like to be able to efficiently generate PSD (positive semi definite) correlation matrices. My method slows down dramatically as I increase the size of matrices to be generated.   1. Could you suggest any efficient solutions? If you are aware of any examples in Matlab, I would be very thankful   2. When generating a PSD correlation matrix how would you pick the parameters to describe matrices to be generated? An average correlation, standard deviation of correlations, eigenvalues?


title: How to conduct conditional Cox regression for matched case-control study?

text: I am attempting to find a program that will let me conduct Cox regression on my matched case-control dataset. Please assist. p.s. I have STATA, SPSS, and MedCalc


title: Self-study question on Confidence Interval

text: Good day everyone, I was making an attempt on a self-study question on Confidence Interval: > The CEO of Micosaft Inc is considering the proposal of offering a child care > program for its employees. As part of the feasibility study, the CEO wishes > to estimate the mean weekly child-care cost of their employees. A sample of > 60 employees who use child care reveals the following statistics: ![enter image description here](http://i.stack.imgur.com/FhkjS.jpg) > where xi is the weekly child-care cost of 1-th employees. Find a 95% > confidence interval for the population mean. **_My attempt as follows:_** ![enter image description here](http://i.stack.imgur.com/VLEgV.jpg) Unfortunately, my answer defers from the model answer of ***80 +_ 2.0111. Appreciate any guidance on where am I not doing right. I've looked at the question and my attempt multiple times without much success to find a way to correct it.


title: Two sample differences between percentage means with negative percentages

text: Hello to everyone and thank you for all the help you've already unknowingly provided to me in the past. I have a group of percentages around zero that are grouped into two categories, and I wanted to know if using a two sample t-test would be appropriate. Both positive and negative values exist, and the values can exceed 100% and -100%, though none of them do by more ~50%. My interest isn't whether the mean is close to zero, but simply if the two groups differ significantly. The data exhibits gaussian distribution, but the mean is not about zero, roughly being -27%. Would there be any problem using a two sample t-test here? I wouldn't assume so, but I haven't worked statistically with percentages, or negative percentages very much, and I just want to make sure.


title: Consequence of negative variance-covaraince matrix on loss

text: I need some help in understanding the impact of a negative varaince-covaraince matrix to loss. For example, I have three loss data sets A,B and C. each has mean and variance. When calculating 3x3 variance-covarance matrix, the determinant is negative. My ultimate goal is to simulate aggregate loss (A+B+C) through assuming distribution etc. I wonder what the consequences of having negative varaince-covaraince matrix to the aggregate loss? Thank you very much!


title: Overly-sensitive goodness of fit measure: where did I go wrong?

text: I'm using the chi-squared goodness-of-fit measure to evaluate instances of a fairly complicated model, but the evaluation is unbelievably sensitive to small variations in the model parameter values. My understanding of the definition and use of this goodness-of-fit measure, and my choice to use it, comes from the book Numerical Recipes, specifically the chapter on data modeling. I work in scientific computing (where this book is a standard reference), not statistics. I'm looking for guidance, as a novice in statistical model evaluation, in what goodness-of-fit measures would be more appropriate to my situation. Here's a simplified example that demonstrates the same sensitivity as I'm seeing in my problem.   * The model is of a "weighted" range on the number line. It has 3 parameters / degrees of freedom: a beginning value, an ending value, and a weight. $(a,b,w)$ is the parameter vector representing the range $[a,b]$ with weight $w$.   * I make 10 observations of any given model instance. Each observation is the weight times the length of overlap between the model instance's range and a particular unit-length interval. My 10 observations are on the intervals $[0,1], [1,2], \ldots, [8,9], [9,10]$. I know a priori that the noise on my observations is Gaussian with unit variance.   * The hidden _ground-truth_ model instance is $(1,9,300)$. So my 10 observations are $$(0.192, 298.839, 301.493, \ldots, 300.273, 300.144, 0.243),$$ which is just $(0, 300, 300, \ldots , 300, 300, 0)$ plus noise. I've got two model instances I want to talk about:   1. The _oracle_ instance $(1,9,300)$ predicts these 10 values: $$(0, 300, 300, \ldots , 300, 300, 0)$$   2. The _candidate_ instance $(1,8.9,300)$ predicts these 10 values: $$(0, 300, 300, \ldots, 300, 270, 0)$$ The chi-squared test statistic, as I'm using it, is: $$\sum_{i=1}^{10} \left(\frac{\mathrm{obs}_i - \mathrm{pred}_i}{\sigma_i} \right)^2$$ where $\sigma_i$ is known to be 1 for all $i$. As I understand it, the chi-squared goodness-of-fit measure for a test statistic $X$ is $$Q(x) = \int_X^\infty f_k(x)\,\mathrm{d}x = 1 - F_k(X)$$ where $f_k$ is the chi-squared PDF with $k$ degrees of freedom, and $F_k$ is the CDF. $k$ is the number of observations minus the number of model degrees of freedom: $k=10-3=7$ in this example. The test statistic for the oracle instance is 9.45, and for the candidate instance it's 918.10. The goodness-of-fit values I get are 0.22 for the oracle (pretty reasonable), and $6{\cdot}10^{-194}$ for the candidate instance. It seems to me that this particular choice of goodness-of-fit measure is way too sensitive: the candidate instance was off by a minuscule value in the parameter space, but the nature of my observations hugely inflates the chi- squared value, and so my goodness-of-fit measure basically ends up screaming: "this model instance is HORRIBLE!" This measure won't provide reasonable values for model instances that are anything but _extremely precisely_ correct. This measure works fine as an objective function for model-fitting: less- correct model instances give even lower goodness-of-fit values. But when I report on the performance of multiple fitting algorithms or models, for example, I want to use the goodness-of-fit values of the best instances I can compute as objective measures by which to compare them. And these values are just laughably small, so tiny that it seems silly to claim that one technique is better than another on the basis of comparing, say, $10^{-194}$ to $10^{-280}$. Their statistical interpretation, also, is a little suspect; this value $Q(x)$ that I compute is a probability of a sort, and probabilities that small just seem meaningless. So: what would be a better choice of measure (or measures) to describe the goodness-of-fit of instances of this model? And, so I can learn from my mistakes: what aspects of the problem should have led me to choose that alternative measure in the first place? [At the risk of making a very long question even longer, here are my answers to questions in the comments:] I have used precisely this value $Q(x)$ (actually, the negative of its log) as my objective function for model fitting. The nature of the problem makes it so that I can only use iterative model refinement to fit, and the example of (1,8.9,300) is analogous to what my fitting process might find after many iterations: not exactly identical to the ground truth, but pretty close, and predicting the vast majority of observations correctly. The essential problem I'm confronting is that, while the oracle instance is definitely a better fit than the candidate instance, the candidate instance was also quite good, and yet the claimed Q(x) is incredibly low. My sense of "quite good" comes from a couple different angles:   * Two of the three parameter values were exactly right, and the third was off by less than 10%.   * Nine of the ten predicted values were correct (well within the bounds, for each observation, of the known observation noise), and the tenth was off by only 10% (though very far down the tail of the noise distribution on the actual observation). It seems like the chi-squared GoF function is overreacting, in a way. I was hoping that this measure would provide some absolute, objective basis for discussing the quality of a fit, but this one seems to reject anything that isn't identical to the ground truth.


title: Replicating bgtest from lmtest

text: _Cross-posted from SO._ I am trying to replicate the results of `bgtest` from the `lmtest` R package. I am using the following dataset:                          rs   month   r20     1    2.365042  1952m3  4.33     2    2.317500  1952m4  4.23     3    2.350833  1952m5  4.36     4    2.451833  1952m6  4.57     5    2.466167  1952m7  4.36     6    2.468417  1952m8  4.11     7    2.485583  1952m9  4.20     8    2.415125 1952m10  4.19     9    2.389875 1952m11  4.15     10   2.418167 1952m12  4.22     11   2.396042  1953m1  4.13     12   2.401042  1953m2  4.10     13   2.400833  1953m3  4.04     14   2.383500  1953m4  3.94     15   2.366708  1953m5  3.95     16   2.365625  1953m6  4.02     17   2.348583  1953m7  3.98     18   2.334375  1953m8  3.94     19   2.133542  1953m9  3.78     20   2.097375 1953m10  3.80     21   2.097708 1953m11  3.78     22   2.130583 1953m12  3.83     23   2.096000  1954m1  3.79     24   2.064042  1954m2  3.79     25   2.115083  1954m3  3.76     26   2.047333  1954m4  3.71     27   1.713875  1954m5  3.65     28   1.606167  1954m6  3.61     29   1.561667  1954m7  3.35     30   1.613292  1954m8  3.36     31   1.621083  1954m9  3.35     32   1.587667 1954m10  3.35     33   1.637792 1954m11  3.38     34   1.865917 1954m12  3.51     35   2.356417  1955m1  3.64     36   3.810000  1955m2  3.85     37   3.797000  1955m3  3.83     38   3.906000  1955m4  4.15     39   3.937000  1955m5  4.21     40   3.969000  1955m6  4.33     41   3.971000  1955m7  4.47     42   4.005000  1955m8  4.84     43   4.072000  1955m9  4.68     44   4.071000 1955m10  4.50     45   4.104000 1955m11  4.64     46   4.072000 1955m12  4.70     47   4.071000  1956m1  4.84     48   5.218000  1956m2  4.87     49   5.165000  1956m3  5.02     50   5.008000  1956m4  4.85     51   4.955000  1956m5  5.12     52   5.136000  1956m6  5.25     53   4.977000  1956m7  5.27     54   4.027000  1956m8  5.20     55   5.091000  1956m9  5.35     56   4.991000 1956m10  5.33     57   5.020000 1956m11  5.50     58   4.858000 1956m12  5.29     59   4.553000  1957m1  4.91     60   4.148000  1957m2  4.93     61   4.099000  1957m3  5.08     62   3.914000  1957m4  5.11     63   3.921000  1957m5  5.43     64   3.854000  1957m6  5.55     65   3.845000  1957m7  5.60     66   4.121000  1957m8  5.75     67   6.605000  1957m9  5.98     68   6.603000 1957m10  5.84     69   6.459000 1957m11  5.89     70   6.375000 1957m12  5.81     71   6.127000  1958m1  5.66     72   6.014000  1958m2  5.65     73   5.523000  1958m3  5.64     74   5.179000  1958m4  5.45     75   4.816000  1958m5  5.46     76   4.294000  1958m6  5.45     77   4.159000  1958m7  5.46     78   3.760000  1958m8  5.49     79   3.625000  1958m9  5.36     80   3.584000 1958m10  5.35     81   3.305000 1958m11  5.36     82   3.152000 1958m12  5.36     83   3.107000  1959m1  5.20     84   3.276000  1959m2  5.20     85   3.287000  1959m3  5.24     86   3.283000  1959m4  5.22     87   3.382000  1959m5  5.28     88   3.452000  1959m6  5.18     89   3.484000  1959m7  5.13     90   3.488000  1959m8  5.21     91   3.472000  1959m9  5.33     92   3.386000 1959m10  5.06     93   3.400000 1959m11  5.04     94   3.687000 1959m12  5.21     95   4.538000  1960m1  5.34     96   4.554000  1960m2  5.43     97   4.621000  1960m3  5.53     98   4.652000  1960m4  5.59     99   4.556000  1960m5  5.62     100  5.681000  1960m6  5.92     101  5.546000  1960m7  5.97     102  5.588000  1960m8  5.95     103  5.565000  1960m9  5.97     104  5.090000 1960m10  5.97     105  4.639000 1960m11  5.97     106  4.349000 1960m12  6.01     107  4.165000  1961m1  6.01     108  4.399000  1961m2  6.04     109  4.485000  1961m3  6.05     110  4.407000  1961m4  6.01     111  4.436000  1961m5  6.08     112  4.537000  1961m6  6.33     113  6.688000  1961m7  6.52     114  6.700000  1961m8  6.63     115  6.552000  1961m9  6.65     116  5.727000 1961m10  6.33     117  5.389000 1961m11  6.34     118  5.403000 1961m12  6.41     119  5.242000  1962m1  6.35     120  5.531000  1962m2  6.26     121  4.405000  1962m3  6.25     122  4.052000  1962m4  6.24     123  3.816000  1962m5  6.25     124  3.921000  1962m6  6.24     125  3.887000  1962m7  5.98     126  3.752000  1962m8  5.77     127  3.635000  1962m9  5.27     128  3.858000 1962m10  5.37     129  3.689000 1962m11  5.42     130  3.717000 1962m12  5.36     131  3.491000  1963m1  5.54     132  3.426000  1963m2  5.74     133  3.756000  1963m3  5.69     134  3.709000  1963m4  5.50     135  3.635000  1963m5  5.31     136  3.702000  1963m6  5.28     137  3.761000  1963m7  5.20     138  3.723000  1963m8  5.22     139  3.674000  1963m9  5.21     140  3.745000 1963m10  5.26     141  3.739000 1963m11  5.51     142  3.721000 1963m12  5.63     143  3.758000  1964m1  5.64     144  4.307000  1964m2  5.85     145  4.302000  1964m3  5.76     146  4.302000  1964m4  5.93     147  4.384000  1964m5  5.90     148  4.464000  1964m6  5.97     149  4.654000  1964m7  6.02     150  4.656000  1964m8  6.00     151  4.703000  1964m9  6.00     152  4.698000 1964m10  6.06     153  6.630000 1964m11  6.23     154  6.627000 1964m12  6.41     155  6.543000  1965m1  6.41     156  6.442000  1965m2  6.43     157  6.549000  1965m3  6.53     158  6.375000  1965m4  6.61     159  6.364000  1965m5  6.76     160  5.542000  1965m6  6.78     161  5.630000  1965m7  6.80     162  5.559000  1965m8  6.65     163  5.559000  1965m9  6.35     164  5.440000 1965m10  6.37     165  5.395000 1965m11  6.40     166  5.521000 1965m12  6.59     167  5.483000  1966m1  6.52     168  5.620000  1966m2  6.61     169  5.604000  1966m3  6.77     170  5.638000  1966m4  6.78     171  5.659000  1966m5  6.82     172  5.728000  1966m6  7.03     173  6.679000  1966m7  7.29     174  6.726000  1966m8  7.41     175  6.747000  1966m9  7.29     176  6.513000 1966m10  6.96     177  6.738000 1966m11  6.97     178  6.527000 1966m12  6.78     179  6.080000  1967m1  6.58     180  6.035000  1967m2  6.49     181  5.495000  1967m3  6.50     182  5.412000  1967m4  6.46     183  5.248000  1967m5  6.65     184  5.275000  1967m6  6.86     185  5.345000  1967m7  6.92     186  5.291000  1967m8  6.90     187  5.475000  1967m9  6.98     188  5.726000 1967m10  7.00     189  7.553000 1967m11  7.22     190  7.484000 1967m12  7.20     191  7.520000  1968m1  7.28     192  7.374000  1968m2  7.28     193  7.108000  1968m3  7.29     194  7.080000  1968m4  7.34     195  7.241000  1968m5  7.50     196  7.242000  1968m6  7.87     197  7.059000  1968m7  7.63     198  6.945000  1968m8  7.63     199  6.577000  1968m9  7.64     200  6.493000 1968m10  7.70     201  6.789000 1968m11  7.93     202  6.777000 1968m12  8.17     203  6.728000  1969m1  8.47     204  7.711000  1969m2  8.61     205  7.782000  1969m3  8.81     206  7.798000  1969m4  8.90     207  7.850000  1969m5  9.46     208  7.880000  1969m6  9.31     209  7.830000  1969m7  9.19     210  7.790000  1969m8  9.49     211  7.811000  1969m9  9.21     212  7.743000 1969m10  8.95     213  7.738000 1969m11  9.29     214  7.650000 1969m12  9.04     215  7.550000  1970m1  9.03     216  7.600000  1970m2  8.79     217  7.270000  1970m3  8.75     218  6.940000  1970m4  8.94     219  6.190000  1970m5  9.40     220  6.870000  1970m6  9.58     221  6.850000  1970m7  9.33     222  6.820000  1970m8  9.19     223  6.820000  1970m9  9.28     224  6.810000 1970m10  9.15     225  6.810000 1970m11  9.51     226  6.820000 1970m12  9.62     227  6.790000  1971m1  9.51     228  6.750000  1971m2  9.35     229  6.660000  1971m3  9.07     230  5.920000  1971m4  9.07     231  5.650000  1971m5  9.03     232  5.590000  1971m6  9.08     233  5.570000  1971m7  9.22     234  5.750000  1971m8  8.96     235  4.830000  1971m9  8.50     236  4.630000 1971m10  8.51     237  4.480000 1971m11  7.79     238  4.360000 1971m12  8.10     239  4.360000  1972m1  7.93     240  4.370000  1972m2  7.90     241  4.340000  1972m3  8.16     242  4.300000  1972m4  8.26     243  4.270000  1972m5  8.60     244  5.210000  1972m6  9.32     245  5.600000  1972m7  9.23     246  5.790000  1972m8  9.36     247  6.440000  1972m9  9.54     248  6.740000 1972m10  9.46     249  6.880000 1972m11  9.45     250  7.760000 1972m12  9.62     251  8.210000  1973m1  9.56     252  8.080000  1973m2  9.65     253  8.070000  1973m3 10.01     254  7.670000  1973m4  9.93     255  7.330000  1973m5 10.02     256  7.060000  1973m6 10.15     257  8.270000  1973m7 10.60     258 10.910000  1973m8 11.30     259 10.970000  1973m9 11.55     260 10.770000 1973m10 11.28     261 11.730000 1973m11 12.00     262 12.460000 1973m12 12.50     263 12.090000  1974m1 12.89     264 11.920000  1974m2 13.50     265 11.950000  1974m3 13.68     266 11.520000  1974m4 14.21     267 11.360000  1974m5 13.80     268 11.230000  1974m6 14.38     269 11.200000  1974m7 14.88     270 11.240000  1974m8 15.29     271 11.060000  1974m9 14.95     272 10.930000 1974m10 15.68     273 10.980000 1974m11 16.75     274 10.990000 1974m12 17.18     275 10.590000  1975m1 16.02     276  9.880000  1975m2 14.58     277  9.500000  1975m3 13.43     278  9.260000  1975m4 13.89     279  9.470000  1975m5 14.53     280  9.430000  1975m6 14.41     281  9.710000  1975m7 13.93     282 10.430000  1975m8 13.87     283 10.360000  1975m9 13.79     284 11.420000 1975m10 14.66     285 11.100000 1975m11 14.81     286 10.820000 1975m12 14.79     287  9.990000  1976m1 13.79     288  8.760000  1976m2 13.46     289  8.460000  1976m3 13.88     290  9.060000  1976m4 13.77     291 10.440000  1976m5 13.59     292 10.960000  1976m6 14.09     293 10.870000  1976m7 14.16     294 10.880000  1976m8 14.33     295 12.050000  1976m9 14.79     296 14.000000 1976m10 16.03     297 14.140000 1976m11 15.79     298 13.780000 1976m12 15.48     299 12.730000  1977m1 14.48     300 11.020000  1977m2 13.93     301  9.920000  1977m3 13.25     302  8.240000  1977m4 13.05     303  7.400000  1977m5 12.69     304  7.450000  1977m6 13.26     305  7.430000  1977m7 13.62     306  6.540000  1977m8 13.12     307  5.680000  1977m9 11.88     308  4.530000 1977m10 10.98     309  4.960000 1977m11 11.28     310  6.370000 1977m12 11.16     311  5.810000  1978m1 11.06     312  5.960000  1978m2 11.75     313  5.930000  1978m3 11.72     314  6.730000  1978m4 12.39     315  8.400000  1978m5 12.72     316  9.170000  1978m6 12.79     317  9.220000  1978m7 12.72     318  8.900000  1978m8 12.55     319  8.980000  1978m9 12.64     320  9.860000 1978m10 12.91     321 11.510000 1978m11 13.16     322 11.570000 1978m12 13.22     323 11.860000  1979m1 13.68     324 12.630000  1979m2 13.94     325 11.350000  1979m3 12.35     326 11.320000  1979m4 11.68     327 11.350000  1979m5 11.94     328 12.570000  1979m6 12.69     329 13.320000  1979m7 12.25     330 13.320000  1979m8 12.30     331 13.380000  1979m9 12.60     332 13.380000 1979m10 13.16     333 15.330000 1979m11 14.54     334 15.900000 1979m12 14.72     335 15.790000  1980m1 14.17     336 16.140000  1980m2 14.45     337 16.180000  1980m3 14.70     338 16.170000  1980m4 14.27     339 16.090000  1980m5 14.01     340 15.800000  1980m6 13.78     341 14.550000  1980m7 13.07     342 14.860000  1980m8 13.58     343 14.400000  1980m9 13.38     344 14.290000 1980m10 13.12     345 13.950000 1980m11 13.22     346 13.070000 1980m12 13.67     347 12.820000  1981m1 13.96     348 12.090000  1981m2 13.89     349 11.530000  1981m3 13.68     350 11.330000  1981m4 13.64     351 11.350000  1981m5 14.31     352 12.090000  1981m6 14.57     353 13.150000  1981m7 15.14     354 13.420000  1981m8 15.09     355 13.960000  1981m9 15.59     356 15.550000 1981m10 15.95     357 14.080000 1981m11 15.44     358 14.510000 1981m12 15.65     359 14.160000  1982m1 15.58     360 13.300000  1982m2 14.74     361 12.480000  1982m3 13.72     362 12.890000  1982m4 13.96     363 12.530000  1982m5 13.69     364 12.230000  1982m6 13.56     365 11.280000  1982m7 13.20     366 10.080000  1982m8 12.23     367  9.910000  1982m9 11.40     368  8.910000 1982m10 10.50     369  9.220000 1982m11 10.64     370  9.960000 1982m12 11.34     371 10.590000  1983m1 11.60     372 10.740000  1983m2 11.50     373 10.470000  1983m3 10.97     374  9.840000  1983m4 10.56     375  9.700000  1983m5 10.65     376  9.470000  1983m6 10.39     377  9.370000  1983m7 10.95     378  9.340000  1983m8 11.07     379  9.160000  1983m9 10.67     380  8.840000 1983m10 10.61     381  8.840000 1983m11 10.29     382  8.870000 1983m12 10.35     383  8.870000  1984m1 10.28     384  8.850000  1984m2 10.42     385  8.430000  1984m3 10.23     386  8.380000  1984m4 10.40     387  8.820000  1984m5 10.93     388  8.860000  1984m6 11.15     389 10.970000  1984m7 11.67     390 10.210000  1984m8 10.98     391 10.020000  1984m9 10.78     392  9.850000 1984m10 10.69     393  9.230000 1984m11 10.32     394  9.100000 1984m12 10.46     395 10.550000  1985m1 10.96     396 12.690000  1985m2 11.06     397 12.930000  1985m3 10.90     398 11.930000  1985m4 10.68     399 11.940000  1985m5 10.88     400 11.890000  1985m6 10.70     401 11.390000  1985m7 10.44     402 10.960000  1985m8 10.37     403 11.060000  1985m9 10.39     404 11.050000 1985m10 10.22     405 11.110000 1985m11 10.37     406 11.150000 1985m12 10.45     407 11.980000  1986m1 10.80     408 12.020000  1986m2 10.40     409 11.060000  1986m3  9.39     410  9.990000  1986m4  8.76     411  9.700000  1986m5  9.00     412  9.320000  1986m6  9.23     413  9.450000  1986m7  9.37     414  9.390000  1986m8  9.41     415  9.610000  1986m9  9.97     416 10.250000 1986m10 10.62     417 10.630000 1986m11 10.80     418 10.660000 1986m12 10.69     419 10.520000  1987m1 10.09     420 10.290000  1987m2  9.83     421  9.350000  1987m3  9.16     422  9.430000  1987m4  9.12     423  8.460000  1987m5  8.82     424  8.540000  1987m6  8.90     425  8.840000  1987m7  9.23     426  9.790000  1987m8  9.20     427  9.690000  1987m9  9.98     428  9.450000 1987m10  9.88     429  8.430000 1987m11  9.20     430  8.190000 1987m12  9.57     431  8.370000  1988m1  9.57     432  8.790000  1988m2  9.38     433  8.270000  1988m3  9.12     434  7.740000  1988m4  9.12     435  7.540000  1988m5  9.27     436  8.880000  1988m6  9.32     437 10.050000  1988m7  9.51     438 11.130000  1988m8  9.47     439 11.530000  1988m9  9.60     440 11.540000 1988m10  9.23     441 12.070000 1988m11  9.30     442 12.540000 1988m12  9.46     443 12.450000  1989m1  9.35     444 12.390000  1989m2  9.15     445 12.410000  1989m3  9.26     446 12.470000  1989m4  9.52     447 12.540000  1989m5  9.52     448 13.590000  1989m6  9.88     449 13.290000  1989m7  9.53     450 13.320000  1989m8  9.37     451 13.440000  1989m9  9.62     452 14.460000 1989m10  9.81     453 14.450000 1989m11  9.99     454 14.500000 1989m12  9.96     455 14.500000  1990m1 10.28     456 14.450000  1990m2 10.72     457 14.570000  1990m3 11.46     458 14.590000  1990m4 11.77     459 14.500000  1990m5 11.49     460 14.380000  1990m6 11.01     461 14.320000  1990m7 11.03     462 14.310000  1990m8 11.41     463 14.260000  1990m9 11.32     464 13.370000 1990m10 11.12     465 12.920000 1990m11 10.94     466 12.960000 1990m12 10.40     467 13.000000  1991m1 10.22     468 12.390000  1991m2  9.89     469 11.640000  1991m3 10.06     470 11.250000  1991m4  9.99     471 10.840000  1991m5 10.15     472 10.720000  1991m6 10.34     473 10.520000  1991m7 10.10     474 10.200000  1991m8  9.89     475  9.660000  1991m9  9.54     476  9.860000 1991m10  9.62     477  9.980000 1991m11  9.68     478 10.100000 1991m12  9.56     479  9.970000  1992m1  9.34     480  9.800000  1992m2  9.21     481 10.100000  1992m3  9.54     482  9.970000  1992m4  9.33     483  9.430000  1992m5  8.99     484  9.420000  1992m6  9.02     485  9.430000  1992m7  8.90     486  9.650000  1992m8  9.13     487  9.160000  1992m9  9.12     488  7.470000 1992m10  9.24     489  6.490000 1992m11  8.84     490  6.390000 1992m12  8.84     491  6.050000  1993m1  8.92     492  5.370000  1993m2  8.63     493  5.380000  1993m3  8.33     494  5.330000  1993m4  8.39     495  5.300000  1993m5  8.60     496  5.190000  1993m6  8.39     497  5.130000  1993m7  7.96     498  5.060000  1993m8  7.39     499  5.170000  1993m9  7.18     500  5.150000 1993m10  7.09     501  4.950000 1993m11  7.06     502  4.870000 1993m12  6.46     503  4.890000  1994m1  6.41     504  4.760000  1994m2  6.83     505  4.830000  1994m3  7.47     506  4.880000  1994m4  7.83     507  4.810000  1994m5  8.24     508  4.880000  1994m6  8.55     509  5.090000  1994m7  8.41     510  5.340000  1994m8  8.52     511  5.390000  1994m9  8.72     512  5.440000 1994m10  8.63     513  5.630000 1994m11  8.53     514  5.870000 1994m12  8.44     515  5.930000  1995m1  8.61     516  6.160000  1995m2  8.52     517  6.090000  1995m3  8.50     518  6.300000  1995m4  8.39     519  6.200000  1995m5  8.18     520  6.370000  1995m6  8.16     521  6.620000  1995m7  8.36     522  6.590000  1995m8  8.24     523  6.520000  1995m9  8.09     524  6.530000 1995m10  8.34     525  6.380000 1995m11  8.01     526  6.220000 1995m12  7.94      which is saved as `ukrates.csv`. Here is the code to attempt to reproduce the `bgtest` module.               rm(list = ls())          library(zoo)     library(lmtest)     library(dynlm)          # read in the data     dfUK = read.csv('./data/ukrates.csv', header = TRUE)     summary(dfUK)          # run the time series regression     zooUK = zoo(dfUK[, c('rs', 'r20')], order.by = as.yearmon(dfUK$month,                                                                    '%Ym%m'))         zooUKAug = merge(zooUK,                           'drs' = diff(zooUK$rs, 1),                       'ldr20' = lag(diff(zooUK$r20, 1), -1))     lmUK2 = dynlm(drs ~ ldr20, data = zooUKAug)          # Breusch-Godfrey regression     zooUKBG = merge(zooUKAug, 'resid' = resid(lmUK2))     lmBG = dynlm(as.formula(paste('resid',                                     '~',                                    attr(lmUK2$terms, 'term.labels'),                                   ' + L(resid, 1)')),                  data = zooUKBG)           # BG test using lmtest package     bgtest(lmUK2, order = 1, type = 'Chisq') # 14.5614          # attempt to recreate BG-test      length(lmBG$residuals)*           sum(lmBG$fitted^2)/sum(lmBG$residuals^2)      This is based on the following code for computing the chi-squared statistic directly from the `bgtest` function code:               > bgtest     function (formula, order = 1, order.by = NULL, type = c("Chisq",          "F"), data = list(), fill = 0)      {         dname <- paste(deparse(substitute(formula)))         if (!inherits(formula, "formula")) {             X <- if (is.matrix(formula$x))                      formula$x             else model.matrix(terms(formula), model.frame(formula))             y <- if (is.vector(formula$y))                      formula$y             else model.response(model.frame(formula))         }         else {             mf <- model.frame(formula, data = data)             y <- model.response(mf)             X <- model.matrix(formula, data = data)         }         if (!is.null(order.by)) {             if (inherits(order.by, "formula")) {                 z <- model.matrix(order.by, data = data)                 z <- as.vector(z[, ncol(z)])             }             else {                 z <- order.by             }             X <- as.matrix(X[order(z), ])             y <- y[order(z)]         }         n <- nrow(X)         k <- ncol(X)         order <- 1:order         m <- length(order)         resi <- lm.fit(X, y)$residuals             Z <- sapply(order, function(x) c(rep(fill, length.out = x),                  resi[1:(n - x)]))             if (any(na <- !complete.cases(Z))) {                 X <- X[!na, , drop = FALSE]                 Z <- Z[!na, , drop = FALSE]                 y <- y[!na]                 resi <- resi[!na]                 n <- nrow(X)             }             auxfit <- lm.fit(cbind(X, Z), resi)             cf <- auxfit$coefficients         vc <- chol2inv(auxfit$qr$qr) * sum(auxfit$residuals^2)/auxfit$df.residual         names(cf) <- colnames(vc) <- rownames(vc) <- c(colnames(X),              paste("lag(resid)", order, sep = "_"))         switch(match.arg(type), Chisq = {             bg <- n * sum(auxfit$fitted^2)/sum(resi^2)                 p.val <- pchisq(bg, m, lower.tail = FALSE)                 df <- m                 names(df) <- "df"             }, F = {                 uresi <- auxfit$residuals             bg <- ((sum(resi^2) - sum(uresi^2))/m)/(sum(uresi^2)/(n -                  k - m))             df <- c(m, n - k - m)             names(df) <- c("df1", "df2")             p.val <- pf(bg, df1 = df[1], df2 = df[2], lower.tail = FALSE)         })         names(bg) <- "LM test"         RVAL <- list(statistic = bg, parameter = df, method = paste("Breusch-Godfrey test for serial correlation of order up to",              max(order)), p.value = p.val, data.name = dname, coefficients = cf,              vcov = vc)         class(RVAL) <- c("bgtest", "htest")         return(RVAL)     }     <environment: namespace:lmtest>      I am wondering why I am getting the different results.


title: How to model time-series data in CRFSuite?

text: I recently came across the CRFSuite package for CRFs. Though, it is primarily used for NLP applications like POS tagging, i was wondering if I could use it to model time-series data as well? Have any of you used CRFSuite? In the tutorials, the raw training data is of the format: He PRP B-NP reckons VBZ B-VP the DT B-NP current JJ I-NP account NN I-NP deficit NN I-NP and they use a file (chunking.py) to convert it to training data format used in CRFSuite. The problem is, a sequence here consists of various labels and CRFSuite learns the model accordingly based on designed features / relationships. However, in my problem, I have a time-series data (force values ) and every value in the sequence is of the same label. For example, I have a time-series of force values of interaction with object type 1 and another time-series of force values of interaction with object type 2, and so inside a sequence all the labels are the same ( corresponding to the object type). One sequence: 1.2 ob1 1.4 ob1 1.5 ob1 1.6 ob1 1.7 ob1 Another sequence: 1.4 ob2 1.3 ob2 1.1 ob2 0.5 ob2 0.1 ob2 I would like to use CRFs to classify objects into object types (ob1 or ob2) by looking at the time-series of force values. If I arrange the training data as shown above, it is unable to capture the relationship (whatever features i decide) and gives very low accuracy on testing data. So, I was wondering if any of you have any insights as to how to represent the data in my domain in this CRFSuite? I have also tried to learn multiple models (each for each category) and then try to classify according to which model gives the highest likelihood of the observed data. But, even that trick did not work. Any pointers? Thanks in advance for your help.


title: Non-Binary Significance Testing?

text: It's been awhile since I've done significance testing... I need to compare 2 average scores: Score A is 7.9 on a 0-10 scale (11pt scale). Score B is 8.2 on a 0-10 scale (11pt scale). The samples sizes are A=70 and B=300. How do I determine if the change in average score is significant or not? Thanks!


title: Good non-informative priors for estimating the parameters of a Gaussian with MCMC (using PyMC)?

text: Say I want to estimate the mean $\mu \in [0, 10] $ of some Gaussian data $\mathbf{x}$ with known variance $\sigma^2 = 1$ using MCMC. Usually I'd use a prior like $\mu \sim \mathrm{Uniform}(0, 10)$ and end up with samples $\hat{\mu}$ that are distributed like this (where $\mu = 5$). This is pretty much what I'd expect, and works fine. Since the samples $\hat{\mu}$ are also Gaussian though, it seems like I should be able to precisely estimate _their_ mean and variance by sampling them directly using priors like these:   * $\mu \sim \mathrm{Norm}(\mu_0, \sigma^2_0)$   * $\sigma^2_0 \sim \mathrm{Uniform}(0, 2)$    * $\mu^2_0 \sim \mathrm{Uniform}(0, 10)$  I hoped to get samples of $\mu_0$ and $\sigma^2_0$ that were very concentrated around the mean ($\approx 5$) and variance ($\approx 0.02$) of $\hat{\mu}$, but instead neither was more precise. $\mu_0$ is spread out a lot and $\sigma^2_0$ is pretty much uniform. Is there a better set of priors/hyperpriors that would be more suitable to sampling in cases like this? I suspect that using uniform, independent hyperpriors might be the issue, but I'm not sure. * * * Here's my PyMC code BTW (possibly there's a bug):               from __future__ import division     import pymc as mc     from pylab import *          # create test data     N = 50     mu = 5     data = randn(N) + mu          # start PyMC variables     mu_0 = mc.Uniform('$\mu_0$', 0, 10)     sigma_0 = mc.Uniform('$\sigma_0^2$', 0, 2)          mu = mc.Normal('$\mu$', mu_0, 1/sigma_0)     data = mc.Normal('data', mu, 1, observed=True, value=data)          # sample     mcmc = mc.MCMC([data, mu, mu_0, sigma_0])      mcmc.sample(iter=50000, burn=5000)          # plot     figure()          for i, v in enumerate(('$\mu$', '$\mu_0$', '$\sigma_0^2$')):         x = mcmc.trace(v)[:].reshape(-1)              subplot(1, 3, i+1)         hist(x, 50)         title(v)          show()


title: Calculating the investment levels for different real estate markets

text: We are trying to figure out how to calculate the real estate investment in different real estate types in different regions. **_We have:_** -The return for the total real estate market. -The return for the individual regional markets (4) (east, south, west, north). -The return for the individual types (5) of real estate (office, apartments, etc.). **_Now we need to calculate:_** -How much of the return in the east/south/west/north is for the individual types of real estate. The reason behind it is that we need to create the best real estate portfolio with regards to the sharpe ratio for example. Regards, Simon Here is our Dataset if you want to have a look at it: http://speedy.sh/2SXNY/Data-Case-2.2.xlsx This paper essentially did exactly what we are trying to do but fails to mention how it is done: http://mitcre.mit.edu/wp- content/uploads/2012/11/HCGPl_REIT_PurePlay.pdf


title: Neural Network hidden layer proportionality

text: I am working on Proportionality on Multiple Hidden layer . My application area is Handwritten Pin code Recognition. I had used five hidden layer for better accuracy. I know time complexity will increase but my work is concentrated on Accuracy only. My network architecture is as follow Input sample are 10 i.e. from 0 to 9 Input layer: 15 rows x 10 column = 150 Neurons Hidden Layer 1 = 102 Hidden Layer 2 = 74 Hidden Layer 3 = 49 Hidden Layer 4 = 37 Hidden Layer 5 = 29 Output layer 16 Neurons because of UNICODE Can anybody help me in pruning this architecture to much extent. I had referred so many papers and afterwards decided to use this network.


title: Why is the restricted boltzmann machine both unsupervised and generative?

text: The restricted boltzmann machine is a generative learning model - but it is also unsupervised? A generative model learns the joint probability `P(X,Y)` then uses Bayes theorem to compute the conditional probability `P(Y|X)`. However, the RBM is an unsupervised feature extractor. There is no `Y`! How can the RBM be said to be an unsupervised algorithm but also generative? Thanks in advance.


title: Detecting parameter influence

text: I have a data set consisting of a system's responses to various test configurations. Every test configuration corresponds to a different parameter set. These parameters can have either continuous values or binary values. I need to figure out which parameters (or which combination) have the most influence upon the system's response. The response is a score between 0-100. What methods can I use for this statistical analysis?


title: How to make linear regression model while 'control' confounding variables?

text: I have data from survey, and Trying to build a linear regression model using R like A~ B however, want to control C, D, E, F, G. like Age, Sex, and other confounding variables. I tried to make some models using lm, glm, etc. and fitted it to my data. However, in Multivariate regression models, I cannot get the graph like below. To use 'segmented' function of R to cut-off point, I guess, I should able to get estimated equation between A and B, while all other variables are controlled. ![enter image description here](http://i.stack.imgur.com/RuWC4.png) The image above describes what I want to do. linear model between A and B, but actual model includes confounders C, D, F, and somehow'controlled' them. The author described that > fitted the two linear regression models for high B and low B, and calculated > the sums of squares of residuals (=observed A -estimated A), from the two > models for each B. The models with the lowest residual sums of squares were > the best models. It seems to be drawn from R(and author said so), but I cannot find any good approach. I'm afraid this question seems to 'tool specific' question/ If so, I'll amend this question, and update it with any help I got. Thanks.


title: Kolmogorov Smirnov conflicts with visual data

text: I recorded the length of 197 nursing care episodes. My normalised data $(x-\mu)/\sigma$ looks like this: ![enter image description here](http://i.stack.imgur.com/gQ4oE.png) Kolmogorov-Smirnov gives me a p-value of 0.14. Can't be right surely?! Am I misinterpreting my results?


title: Residual from Functional Demographic model for fertility

text: I'm trying to fit a functional demographic model (`fdm`) to fertility rates using the `demography` package in R. When I've plotted the residuals using `filled.contour`, it seems that residuals persist in the lowest age interval over the period. Since residuals should normally be random, is it normal that the residuals are behaving like this in this case or does it mean my model is underestimating fertility rate in this age group? Or am I simply misinterpreting the graph? Note that I've also tried the `adelaid.R` codes for Australian fertility given by Prof. Rob Hyndman, and it's same outcome. Here is the code:               library(demography)         library(addb)      ausfert <- extract.years(aus.fert, years = 1921:2000)       ausf.sm <- smooth.demogdata(ausfert)       ausf.fit <- fdm(ausf.sm,order = 3, method = "classical")       plot(residuals(ausf.fit), type = "filled.contour")      ![The codes resulted in this graph](http://i.stack.imgur.com/QShny.png)


title: r: normalmixEM not consistent?

text: I have the following data: http://s000.tinyupload.com/?file_id=00083355432555420222 I want to fit a mixture density of two normal distributions. I have the formula: \begin{align} f(l)=\pi \phi(l;\mu_1,\sigma^2_1)+(1-\pi)\phi(l;\mu_2,\sigma^2_2) \end{align} my R code is:               normalmix<-normalmixEM(dat,k=2,fast=TRUE)          pi<-normalmix$lambda[1]         mu1<-normalmix$mu[1]     mu2<-normalmix$mu[2]         sigma1<-normalmix$sigma[1]     sigma2<-normalmix$sigma[2]      Now I have the problem, that the output is not consistent, i.e. every time I run the code, I get different outputs! And they are very different, no small differences, which could be due to the precision of the numerical procedures. E.g. sometimes for pi I get               [1] 0.2653939      or               [1] 0.3318069      I already recognized, that sometimes the numbering is changed, so the pi of 0.7 would be equal to a pi of 0.3. Okay, I got this, I don't know why the R procedures does this, but this would not be a problem. But the problem is, that the outputs are way to different, sometimes I even get an error message (german): Fehler in while (dl > eps && iter < maxit) { : Fehlender Wert, wo TRUE/FALSE nötig ist Also, the number of iterations is very different, from 29 up to 1000 ...... edit: I now tried to set the initial values manually.... In the Rmanual it says: > Lambda: "Initial value of mixing proportions. Automatically repeated as > necessary to pro- duce a vector of length k , then normalized to sum to 1" I started to try it with giving different values of lambda: 1.               normalmix<-normalmixEM(dat,lambda=c(0.5,0.5),k=2,fast=TRUE)     normalmix$loglik         normalmix$mu     normalmix$sigma      which gives a loglik of `617.2996` and the parameters:               > normalmix$mu         [1]  0.0003769442 -0.0008892282         > normalmix$sigma     [1] 0.02814997 0.01232638      2.               normalmix<-normalmixEM(dat,lambda=c(0.8,0.2),k=2,fast=TRUE)     normalmix$loglik         normalmix$mu     normalmix$sigma      which also gives a loglik of `617.2996` ? and the parameters               > normalmix$mu         [1]  0.0003767148 -0.0008892515         > normalmix$sigma     [1] 0.02814819 0.01232488      mh, I don' know.... then, I tried it with different values of $\mu$: 1.               normalmix<-normalmixEM(dat,mu=c(-0.001,0.002),k=2,fast=TRUE)     normalmix$loglik         normalmix$mu     normalmix$sigma      which also leads to a loglik of `617.2996` and the parameter output:               > normalmix$mu         [1]  0.0003767138 -0.0008892516         > normalmix$sigma     [1] 0.02814818 0.01232487      2.               normalmix<-normalmixEM(dat,mu=c(0.02,0.01),k=2,fast=TRUE)     normalmix$loglik         normalmix$mu     normalmix$sigma      which again! leads to a loglik of `617.2996` and the parameter output:               > normalmix$mu         [1]  0.0003769462 -0.0008892280         > normalmix$sigma     [1] 0.02814999 0.01232639     >       what should I do now?


title: What is the “Logistic Regression”? I cannot have a unified concept

text: I'm interested in logistic regression for modeling a classification problem. I tried to study logistic regression with two books, "Discrete Choice Methods with Simulation" (Train) and "Applied Logistic Regression" (Hosmer and Lemeshow). To clarify my question, I picked two logit model examples from those books.   1. > Household's choice between a gas and an electric heating system (binary logit)   >  Let the subscripts $g$ and $e$ denote gas and electric, $PP$ and $OC$ are > the purchase price and operating cost, $\beta_1$ and $\beta_2$ are scalar > parameters. Then the probability that the household choose gas heating is   >  $$P_g = > \frac{e^{\beta_1PP_g+\beta_2OC_g}}{e^{\beta_1PP_g+\beta_2OC_g}+e^{\beta_1PP_e+\beta_2OC_e}}$$   2. > Relationship between age and the presence or absence of CHD (binary logit)   >  Let $x$ be the age, the subscripts $0$ and $1$ denote absence and presence > of the CHD and $\beta_1^0$, $\beta_2^0$, $\beta_1^1$ and $\beta_2^1$ are > scalar parameters. Then the probability that a $x$ years old person has a > CHD is:   >  $$P_1 = > \frac{e^{\beta_1^1+\beta_2^1x}}{e^{\beta_1^1+\beta_2^1x}+e^{\beta_1^0+\beta_2^0x}}$$   >  (Although the model needs two scalar parameters to be estimated, I wrote > the formula in this form to clarify the question.) I am familiar with the logistic regression showed in the second example. The independent variable (age in this example) is included in both choices (or classes), and of course, the value of independent variable is same for different classes. And there are coefficients for each class and each independent variable which has to be estimated, and those coefficients are different for different classes. However, the model of the first example has variables which are different for each choices ($PP_e$, $PP_g$, $OC_e$ and $OC_g$), and the coefficients are the same between choices. Roughly, I can understand two models separately, but I cannot build the unified concept of the logistic regression. They seems different to me. For example, many textbooks and materials provide the (multinomial) logistic regression formula, which contains independent variables and the coefficients (which are different for each classes). With this formula, I can understand the second example, but can this formula incorporate the situation showed in the first example? In other words, does this formula have the variables which are different for each class?


title: Statistical test to investigate multiple data sets whether they are independent

text: I had an 18 variable original data. I did principal component analysis and converted it to 6 PCs explaining 75% of the variance of original data. The I clustered the 6 PCs. For that, I first used average linkage to calculate the mean value of the PC scores with certain number of clusters (I tried 3,4,5,6). Then I used the mean values and number of clusters as seeds for K means clustering and got 4 clusters. Now I have got a problem. What test should i use to know how significant is the clustering procedure? one way is of course calculating the means and standard deviation of the original 18 variables with each of the four clusters. But I am looking for a statistical test which can prove that the new 4 data sets are independent and the clustering is significant.


title: correlation coefficient

text: Hi I am finding a correlation between a count variable (var1) (number of errors) and an ordinal variable (var2) (never, weekly, monthly, yearly). I have tried the following commands in Stata:               twoway scatter var1 var2     spearman var1 var2, stats(rho obs p)      Not sure what if I am going in the right direction? Please can some one have any thoughts!


title: Best infographic of statistics as a subject

text: I am trying to figure out the complete ecosystem of statistics theory. **A Classification Tree of all the math/stat theories over time** Has anybody done or have any links to Infographics about statistics. perhaps starting with stats based on past, present or future. eg, billing, current price of stock, prediction of future earnings... or different styles, induction, hypothesis testing, Bayesian... Maybe Cluster analysis by industry or types of datasets. eg. Marketing and Advertising, Basic time series for measuring clicks Optimisation or Principal Component Analysis for creative optimization, I don't know, something where somebody new can figure out the scope of statistics. I mean I know it can be used in pretty much every walk of life, but there must be some sort of general classifications. Something like this but for statistics. http://static.bbci.co.uk/naturelibrary/3.1.5//images/tree_of_life_full.gif Please feel free to edit this post if you can improve the question. Thanks


title: Three important issues

text: I don't understand three important issues : (1) Why do we calculate the test statistic under the null hypothesis ? (2) Why do we reject the null hypothesis when : p-value < significance level ? (3) Why is type-I error more serious than type-II error ?


title: Prior for Bayesian Inference on Failure Rate in Poisson Distribution

text: I'm trying to derive the posterior distribution for the failure rate (lambda) of a process with poisson distribution. I have tried the use of an improper uniform distribution on lambda by letting the max tend to infinity. My posterior then just has the Poisson like form again (but as a likelihood function of continuous variable lambda). I've also tried using various expert opinions as confidence intervals - ranging from what I consider optimistic to very conservative - and approximating by a gamma distribution. From these approaches I get very similar results, but I'm not sure I have a sufficiently good argument to defend my use of priors. So, I'm concerned my approach is not robust and on reading about Jeffrey's prior I'm wondering if I should learn more about it and how to use it. Any advice?


title: I'm getting "jumpy" loadings in rollapply PCA in R. Can I fix it?

text: I have 10 years of daily returns data for 28 different currencies. I wish to extract the first principal component, but rather than operate PCA on the whole 10 years, I want to rollapply a 2 year window, because the currencies' behaviours evolve and so I wish to reflect this. However I have a major problem, that is that both the princomp() and prcomp() functions will often jump from positive to negative loadings in adjacent PCA analyses (ie 1 day apart). Have a look at the loading chart for the EUR currency: ![enter image description here](http://i.stack.imgur.com/8A2oH.png) Clearly I can't use this because adjacent loadings will jump from positive to negative, so my series which uses them will be erroneous. Now take a look at the absolute value of the EUR currency loading: ![enter image description here](http://i.stack.imgur.com/KMkhE.png) The problem is of course that I still cannot use this because you can see from the top chart that the loading does go from negative to positive and back at times, a characteristic which I need to preserve. Is there any way I can get around this problem? Can I force the eigenvector orientation to always be the same in adjacent PCAs? By the way this problem also occurs with the FactoMineR PCA() function. The code for the rollapply is here:               rollapply(retmat, windowl, function(x) summary(princomp(x))$loadings[, 1], by.column = FALSE, align = "right") -> princomproll


title: Bayesian Linear Regression is so hard to understand?

text: I'm learning Bayesian Linear Regression from a book, the linear model is $$p(w|x,\phi,\sigma^2)=Norm_w[\phi^Tx,\sigma^2]$$, as put in the book, we use Bayes approach to do the parameters estimation. Here comes the problem: I thought it is pretty clear that we should introduce a **conjugate** prior for the parameters $\phi$ and $\sigma^2$, which should have a **normal-scaled inverse gamma distribution** , right? But the book first assume that $\sigma^2$ is known, and introduce a prior distribution for $\phi$ alone, which is a 0-mean Gaussian, and does the estimation for $\phi$. After all that, it assumes $\sigma^2$ is not known, and estimates it. Why should we separate them? UPDATE The book is _Machine Learning: A Probabilistic Perspective_ , p. 232, Section "Baysian Linear Regression." I just found this article which also assumes one is known and later assumes it is unknown. Bayesian Linear Regression


title: How to do dimensionality reduction in R

text: I have a matrix where a(i,j) tells me how many times individual i viewed page j. There are 27K individuals and 95K pages. I would like to have a handful of "dimensions" or "aspects" in the space of pages which would correspond to sets of pages which are often viewed together. My ultimate goal is to then be able to compute how often individual i has viewed pages that fall in dimension 1, dimension 2, etc. I have read the R documentation on principal component analysis and single value decomposition and have executed these commands, but I am unsure how to proceed. How can I use dimensionality reduction to do this? Or is this really a clustering problem and I should instead look into clustering algorithms? Many thanks for any insight ~l


title: Does this single value match that distribution?

text: this feels like a very naive question but I'm having difficulty seeing the answer. I have one set of 30 values. Independently I obtained a 31st value. Null hypothesis is that the 31st value is part of the same distribution. Alternative is that its different. I want some kind of p-value or likelihood measure. Some thoughts I've had:   * This is similar to wanting to do a two-sample t-test - except that for the second sample I only have a single value, and the 30 values aren't necessarily normally distributed.   * If instead of 30 measurement I had 10000 measurement, the rank of the single measurement could provide some useful information.  How can I calculate this likelihood or p-value? Thanks! Yannick


title: Why are after training 5-state HMM only few entries of transition matrix left greater than zero?

text: I try to create the speech recognition system based on 5-state HMM + Multivariate Gaussian function. I use my own feature vector derived from MFCC (Mel-frequency cepstral coefficients). The problem is that after training only a few entries of the HMM transition matrix are left greater then zero. It should be mentioned that only states 2 -- 4 can emit observation. States 1 and 5 help connecting model with other models to make a chain. The initial matrix: $$\left[ \begin{array}{ccccc} 0.0 & 1.0 & 0.0 & 0.0 & 0.0 \\\ 0.0 & 0.6 & 0.4 & 0.0 & 0.0 \\\ 0.0 & 0.0 & 0.6 & 0.4 & 0.0 \\\ 0.0 & 0.0 & 0.0 & 0.7 & 0.3 \\\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \end{array} \right]$$ Matrix after model training on phoneme "a" samples is looking like that: $$\left[ \begin{array}{ccccc} 0.0 & 1.0 & 0.0 & 0.0 & 0.0 \\\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0 \\\ 0.0 & 0.0 & 0.8910256 & 0.1089744 & 0.0 \\\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0 \\\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \end{array} \right]$$ The training corpus consists of 187 sequences with about 10 vectors in average in every sequence. The performance of such model is quite poor. So what might be the main reasons for such sparse transition matrix after model training? Thank you for any help in advance.