title: How to make the rare events corrections described in King and Zheng (2001)?

text: I have a dataset with a binary (survival) response variable and 3 explanatory variables (`A` = 3 levels, `B` = 3 levels, `C` = 6 levels). In this dataset, the data is well balanced, with 100 individuals per `ABC` category. I already studied the effect of these `A`, `B`, and `C` variables with this dataset; their effects are significant. I have a subset. In each `ABC` category, 25 of the 100 individuals, of which approximately half are alive and half are dead (when less than 12 are alive or dead, the number was completed with the other category), were further investigated for a 4th variable (`D`). I see three problems here:   1. I need to weight the data the rare events corrections described in King and Zheng (2001) to take into account the approximate 50% - 50% is not equal to 0/1 proportion in the bigger sample.    2. This non-random sampling of 0 and 1 leads to a different probability for individuals to be sampled in each of the `ABC` categories, so I think I have to use true proportions from each category rather than the global proportion of 0/1 in the big sample.   3. This 4th variable has 4 levels, and the data are really not balanced in these 4 levels (90% of the data is within 1 of these levels, say level `D2`). I have read the King and Zheng (2001) paper carefully, as well as this CV question that led me to King and Zheng (2001) paper, and later this other one that led me to try the `logistf` package (I use R). I tried to apply what I understood from King and Zheng (2001), but I am not sure what I did is right. I understood there are two methods:   * For the prior correction method, I understood you only correct the intercept. In my case, the intercept is the `A1B1C1` category, and in this category survival is 100%, so survival in the big dataset and the subset are the same, and therefore the correction changes nothing. I suspect this method should not apply to me anyway, because I do not have an overall true proportion, but a proportion for each category, and this method ignores that.   * For the weighting method: I calculated _w i_, and from what I understood in the paper: "All researchers need to do is to calculate _w i_ in Eq. (8), choose it as the weight in their computer program, and then run a logit model". So I first ran my `glm` as:              glm(R~ A+B+C+D, weights=wi, data=subdata, family=binomial)      I am not sure I should include `A`, `B`, and `C` as explanatory variables, since I normally expect them to have no effect on survival in this subsample (each category contains about the 50% dead and alive). Anyway, it should not change the output a lot if they are not significant. With this correction, I get a good fit for level `D2` (the level with most of individuals), but not at all for others levels of `D` (`D2` preponderates). See the top right graph: > ![Fits](http://i.stack.imgur.com/AO3ml.png) Fits of a non-weighted `glm` > model and of a `glm` model weighted with _w i_. Each dot represents one > category. `Proportion in the big dataset` is the true proportion of 1 in the > `ABC` category in the big dataset, `Proportion in the sub dataset` is the > true proportion of 1 in the `ABC` category in the subdataset, and `Model > predictions` are the predictions of `glm` models fitted with the subdataset. > Each `pch` symbol represents a given level of `D`. Triangles are level `D2`. Only later when seeing there is a `logistf`, I though this is perhaps not that simple. I am not sure now. When doing `logistf(R~ A+B+C+D, weights=wi, data=subdata, family=binomial)`, I get estimates, but the predict function does not work, and the default model test returns infinite chi squared values (except one) and all p-values = 0 (except 1). Questions:   * Did I properly understand King and Zheng (2001)? (How far am I from understanding it?)   * In my `glm` fits, `A`, `B`, and `C` have significant effects. All this means is that I deparse a lot from the half / half proportions of 0 and 1 in my subset and differently in the different `ABC` categories – isn't that right?   * Can I apply King and Zheng's (2001) weighting correction despite the fact that I have a value of tau and a value of $\bar y$ for each `ABC` category instead of global values?   * Is it an issue that my `D` variable is so unbalanced, and if it is, how can I handle it? (Taking into account I will already have to weight for the rare event correction...Is "double weighting", i.e. weighting the weights, possible?) Thanks! **Edit** : See what happens if I remove A, B and C from the models. I do not understand why there is such differences. > ![Fits2](http://i.stack.imgur.com/gCySf.png) Fits without A, B, and C as > explanatory variables in models


title: Do unique visitors to a website follow a power law?

text: Suppose I have an ordered vector where the first element is the number of visits to a website in a given period of time by the unique IP with the highest number of visits, the second element is the number of visits by the unique IP with the second highest number of visits, and so on. I understand there may be per site variations, but is there in general an assumed pattern to the shape of this vector? Does it, for instance, follow a power-law distribution?


title: How to calculate error matrices, K hat and var(K hat)?

text: I have a very simple data set consisting of three columns: Ground Truth Canopy Class, Method 1 Canopy Class and Method 2 Canopy Class. Each row in the columns represents the canopy class (i.e. 1 through 5). I have produced an error matrix in excel and calculated the overall accuracy and Khat (Figure 1). However, now I need to test whether or not method one is significantly different from method two using the attached equations (Figure 2 and 3). I could use help to calculate the following using R or Excel:   1. Z-score    2. K hat   3. Variance of K hat Is there a package in R I've overlooked that can aid me these calculations? Figure 1 ![enter image description here](http://i.stack.imgur.com/mQNIy.jpg) Figure 2 ![](http://i.stack.imgur.com/hRdJs.jpg) Figure 3 ![](http://i.stack.imgur.com/U0r9u.jpg) Source: Assessing the Accuracy of Remotely Sensed Data (Congalton, 2009)


title: Is using an ANOVA or multiple t-tests the best way to compare monthly and annual data?

text: I have daily precipitation data for two locations over two years. I need to compare them to see if there is an annual difference (combining the two years of data), and if there is a difference in any particular month (a month is also the combination of the two years of data). Is this best analyzed using a 2 x 12 two-way ANOVA (although N is different for different months), or would a lot of student t-tests be what I need, or both?


title: Is there an alternative to the Kolmogorov-Smirnov test for tied data with correction?

text: I've got a bunch of data from two samples (control and treated), each containing several thousand values which are to undergo significance testing in R. Theoretically, the values should be continuous, but due to rounding done by measurement software they aren't and they have got ties. The distributions are unknown and the shapes of control and treated distributions might be different, so I'd like to use a non-parametric test to compare if the difference across the samples is significant for 10 different factors. I thought of using the Kolmogorov-Smirnov test, but it's not really suitable for ties. I recently stumbled upon a new R library called Matching that executes a bootstrap version of K-S test and tolerates ties. Now is this really a good idea or should I use another test instead? And do I need to adjust the p-value?


title: I got a good OLS fit for integer variables, do I still need to use count data methods?

text: First of all, I'm not a statistician, but I'm teaching myself some methods I require for a project I'm doing now. I have a 2D dataset of N observations. For the ith observation, the first entry is the estimated number of lines in the ith chapter of a book, and the second entry is the actual number of lines in that chapter. All the values are large (so I guess there's no point on using Poisson distributions), and after taking log to both entries I get a nice line. There is no heteroscedasticity, and the distribution of residuals seems to be normal. My question is, should I anyway necessarily use methods for count data? If not, how could I justify it? If yes, a negative binomial distribution would do the trick?


title: Newey-West t-statistics

text: I have a time-series which is autocorrelated by construction, and might be heteroscedastic. I have calculated the sample mean of this time-series, and would like to calculate the t-statistic corresponding to the hypothesis that the mean of this time-series is zero. It is my understanding that since my time-series is autocorrelated and possibly heteroscedastic, that I must use a t-statistic "adjusted for serial dependence according to the Newey-West method".   1. I have problem understanding the method, and how to implement this in Matlab. As far as I understand, Newey-West is used in regressions to obtain HAC standard errors, since the OLS standard errors are not a reliable basis for inference under serial correlation of the error term in a regression. But in my case, I am not regressing anything, so how does the Newey-West method fit in? Under no autocorrelation and homoscedasticity, I would have simply divided the sample mean (minus 0) by the standard error (sample standard deviation over square root of the number of observations). There would therefore be no need for regressing anything.    2. How can I implement computation of Newey-West t-statistics in Matlab?


title: Google trends data for interest

text: I was discussing about the popularity of some terms and used google trends to conclude in the decrease of their popularity. Here is an exemple of the queries for some of the biggest french engineering schools. I concluded rapidly in a decreasing interest/popularity of these schools. But giving another thought, I wondered if the decrease may have another explanation. Google trends give a relative number, proportional to the total number of queries. Is it possible that the population accessing google has so much evolved in these years that we can't conclude to a decrease in interest ? I mean that the majority of queries should come from educated/wealthy people and researchers. So with the democratization of internet access the fraction of wealthy/educated people over total should have decreased. Thus explaining a decrease of the querries over total. How can we show this effect ? How can we remove this effect ?


title: Now that I've rejected the null hypothesis what's next?

text: I've time and again _rejected_ or _failed to reject_ the null hypothesis. In the failure to reject case, you conclude that there isn't sufficient evidence for rejection and you "move on" (i.e., either gather more data, end the experiment etc.,) But when you "do" reject the null hypothesis, providing _some_ evidence for the alternative hypothesis you can't really "prove" that your alternative hypothesis indeed holds true. So, what are the common next steps once you reject the null hypothesis? What tools/techniques does one adopt to "analyze the problem further" to make more the findings more conclusive? What are the logical "next steps" as a statistician warranting further analysis? For example: $H_0: \mu_1 = \mu_0$ $H_1: \mu_1 > \mu_0$ (say we know the expected direction) Once we reject the null hypothesis at some level of significance we have "some evidence" for the alternative to be true, but we can't draw that conclusion. If I really want to draw that conclusion conclusively (pardon the double word play) what should I do? I've never pondered this question during my undergrad days but now that I'm doing a fair deal of hypotheses testing I can't help but wonder what's ahead :)


title: Analytic or sample standard deviation with binomial data

text: I've been looking for recommendations on whether it's better to use the sample standard deviation (SD) for a binomial distribution or use the analytic SD (or variance). It's for experiments with accuracy data so each subject contributes several binomial (correct or incorrect) responses resulting in an accuracy score between 0 an 1 (Ncorrect / Ntotal). I get the mean across subjects. Assuming each one contributes the same amount should I used the sample SD of the individual accuracy scores or use the analytic equation (`sqrt(pq / n`) where `n` is the number of responses / subject)? This is a query for descriptive purposes only. I could use a multi-level logistic regression to model but I'm just looking for simple descriptives. **Example:** Each subject gets 20 tries at the task and there are 6 subjects with accuracy `0.95, 0.80, 0.80, 0.65, 0.90, 0.70`. The mean accuracy would be `0.8` and therefore the analytic SD of that accuracy is `sqrt(0.8 * 0.2 / 20) = 0.089`. However, the SD of those six numbers calculated on the sample is 0.114. Which is the better SD estimate to use?


title: Is Rule of Three inappropriate in some cases?

text: I'm building binomial proportion confidence intervals for a patient dataset containing the frequency of home nursing visits during the week prior to hospital admission. The freq. of home visit categories are integer values of 1, 2, 3,... and observations for a particular patient can be assumed to be independent from other patients. It looks something like this for a dataset of 100 patients:               No. Visits:   Frequency:      % Total:         1             10             10%         2             50             50%         3             30             30%         5             10             10%     Total:           100            100%       So in this example I'm building separate 95% CIs for $p_1 = 0.10$, $p_2 = 0.50$, etc. with $n=100$. I'm using the Rule of Three to compute 95% CI $= [0, 3/n]$ for frequency of visits (in this case No. Visits $= 4$) which were not observed in the dataset. However it occurred to me I could mechanically use the Rule of Three to compute 95% CIs $= [0, 3/n]$ extending for No. Visits $= 6$ all the way up to infinity. Is this reasonable?


title: Error Correction Model in Stata

text: I have run an ECM in Stata, where my model has six variables, having found 3 cointegrating relationships.   When putting this table into a paper, what are the coefficients I need to report on each variable and where is the error correction term? Are there multiple error correction terms given that there are 3 cointegrating relationships?


title: Which coefficient: correlation or regression?

text: How can I interpret the results if the inter-item correlation is negative while regression coefficient is positive? I have run correlation to ensure if I can work on regression. Where should be my point of discussion be focused? On the correlation or regression?


title: Data has two trends; how to extract independent trendlines?

text: I have a set of data that is not ordered in any particular way but when plotted clearly has two distinct trends. A simple linear regression would not really be adequate here because of the clear distinction between the two series. Is there a simple way to get the two independent linear trendlines? For the record I'm using Python and I am reasonably comfortable with programming and data analysis , including machine learning but am willing to jump over to R if absolutely necessary. ![enter image description here](http://i.stack.imgur.com/SkEm3.png)


title: Assumptions of linear models and what to do if the residuals are not normally distributed

text: I am a little bit confused on what the assumptions of linear regression are. So far I checked whether:   * all of the explanatory variables correlated linearly with the response variable. (This was the case)   * there was any collinearity among the explanatory variables. (there was little collinearity).   * the Cook's distances of the datapoints of my model are below 1 (this is the case, all distances are below 0.4, so no influence points).   * the residuals are normally distributed. (this may not be the case) But I then read the following: > violations of normality often arise either because (a) the distributions of > the dependent and/or independent variables are themselves significantly non- > normal, and/or (b) the linearity assumption is violated. **Question 1** This makes it sound as if the independent and depend variables need to be normally distributed, but as far as I know this is not the case. My dependent variable as well as one of my independent variables are not normally distributed. Should they be? **Question 2** My QQnormal plot of the residuals look like this: ![normality check of residuals](http://s22.postimg.org/b1wevtw8h/p1_normality_residuals.png) That slightly differs from a normal distribution and the `shapiro.test` also rejects the null hypothesis that the residuals are from a normal distribution:               > shapiro.test(residuals(lmresult))     W = 0.9171, p-value = 3.618e-06      The residuals vs fitted values look like: ![residuals vs fitted](http://s28.postimg.org/4o737gi3x/p1_residualsvsfitted.png) What can I do if my residuals are not normally distributed? Does it mean the linear model is entirely useless?


title: An R function for performing searches

text: I am slightly new to writing functions in R. Here I have a basic function that searches for a pattern and returns the indexes of where the occurence occurs given a list dataset               #this functions takes a pattern and prints the indexes for the matches     find_domain <- function(pattern,list) grep(pattern,list,ignore.case=T)      Sometimes during the analysis i found that i need to do that for a larger dataset. The idea is to look for occurence of a domain (pattern) in each of the d7_dataset rows(since each row represents a protein sequence and each column is a domain) such that i can then count the number of sequences that contain or have a given domain(pattern in this case) So i wrote this               seqs <- c()     for (i in 1:nrow(d7_dataset)){            pos <- find_domain("CIDRγ13",d7_dataset[i,])            if (!is.na(pos[1]))              seqs <- c(seqs,1)            }          total_seqs <- sum(seqs)     total_seqs       However this seems like it can be written as a single function so that i can eliminate duplication and simplify it such that it is generic enough to apply to multiple datasets. Any ideas on how to condence it to a single function? In a nutshell, Given a pattern or list of patterns and a dataset(as a dataframe) to search, return the number of sequences(rows in the dataset) that contain the given pattern. More or less like this               results <– search(a list of patterns to search,dataset)      results should look like               pattern  count      pat1     40      pat2      3      pat3      0      .         .       .         .      It would be very good if the pattern argument would accept a list to search Thank you


title: how to calculate correlation when a variable number of observations are collected for each subject?

text: I have collected observations (y variable) from people of various ages (x variable). Some people provided multiple data points (anywhere from 2-10), whereas others provided just one. I want to see if x is correlated with y. How would I account for the fact that some people gave more observations than others? This is not a time series dataset. Thanks for any help!


title: How to calculate likelihood for a bayesian model?

text: I am trying to do Bayesian posterior predictive checking, whereby I calculate the DIC for my fitted model, and compare to DIC from data simulated from the fitted model. I can get the DIC out of winBUGS, however I am not sure how to calculate the likelihood (for the DIC) outside of winBUGS (i.e., without fitting new models). All of the literature is quite general with respect to discrepancy functions for model checking, and have notation like p(D|theta) which I understand, but doesn't help me when I want to actually do the calculation. Given a set of parameter values, and a Y_rep data set - how do I calculate the likelihood?


title: LMG don't sum to model's R2

text: I have a regression model and I want to determine the relative importance of the predictors. I used the package `relaimpo` in r. The package says that 98.84% of the variance is explained by the model               Response variable: C      Total response variance: 115.4857      Analysis based on 161187 observations           10 Regressors:      I(VPT^(-0.66)) TPA S SD TPA:VPT S:SD TPA:SD TPA:S SD:VPT S:VPT      Proportion of variance explained by model: 98.84%      I did not normalize the metrics.               Metrics are not normalized (rela=FALSE).       These are my relative importance metrics:               Relative importance metrics:                                   lmg     I(VPT^(-0.66)) 0.2647988336     TPA            0.0776291924     S              0.0816760485     SD             0.2661375234     TPA:VPT        0.0779008229     S:SD           0.0036628215     TPA:SD         0.0002298444     TPA:S          0.0002684624     SD:VPT         0.0065747570     S:VPT          0.0114800491      They sum up `0.79`. Shouldn't they sum up to `0.9884`? If so, why don't they?


title: Digging deeper into the FitDist function

text: I am using the script `allfitdist` to find which distribution best fits my data according to the tests included (BIC, AIK etc...) - these all stem from the script using the `fitdist` function in the statistics toolbox. I manage to do this and get a nice historgram with the various distributions superimposed on my data. This question poses something similar to what I am struggling with, but I want to go a step further and know if it is possible to get the further parameters (for example the skew and kurtosis) of the output distributions. I have literally no idea how I can do this. The output for a norm dist only contains mean and st.dev, but it seems from the figure that there is skewness and kurtosis - I wouldn't believe them to both be exactly zero either, given my data is market return of a stock (unless the `fitdist` function used has that as a assumption??) If it would be helpful, I could post a picture of the output.


title: How to estimate a critical value of Spearman's correlation for n=100?

text: I have searched a lot, and I can only find tables that show critical values up to n=30. Can someone provide, or point me to, a simple method of estimating this value for different $\alpha$?


title: Kolmogorov-Smirnov and lattice paths

text: I have seen asserted that the problem of computing the null distribution of Kolmogorov's $D_n^+$ statistic for a finite sample size maps onto the problem of computing the number of lattice paths that stay below the diagonal, and thus can be solved by the ballot theorem. I am familiar with lattice paths and the ballot problem. I am also familiar the expression of the distribution of $D_n^+$ as a series of integrals. But I don't see how one problem maps onto the other. Can someone explain or point me to an ariticle or book that does? I also see the claim the the null distribution of the Kolmogorov-Smirnov $D_n = \max(D_n^+,D_n^-)$ maps onto another lattice path problem that could be solved by a "two-sided ballot theorem". I don't know what a "two-sided" version of the ballot problem would be. Again, can someone explain or point me to an explanation? Finally, is there a general framework around all of this? Can the Kuiper statistic be mapped to yet another lattice path problem? The two-sample KS test? The AD statistic?


title: Two-stage clustering in R

text: * Is it possible to do 2-stage cluster analysis in R?   * Can anybody provide me resource on it?


title: Reality check using GLM

text: I need a reality check if you will. I have a data set, where I know how many individual butterflies of two species co-occur at one meadow (not always, though). I have additional variables, for instance wet/dry meadow, intensely cultivated/not cultivated, percent of area around the meadow covered by wet or dry meadows... This is the head of the dataset. 50 rows in total, 25 per species. Notice that all columns are identical except count and species, indicating that they come from the same sampling location.               > head(dej)       count     type1 type2 perc.for.100m perc.dry.100m perc.wet.100m species     1     1 intensive   dry        13.836        22.724         0.000   reali     2     3 extensive   wet         6.877         1.613        52.213   reali     3     4 intensive   wet        22.770         0.537        44.901   reali     4     6 intensive   dry        17.346        42.322         6.359   reali     5     1 extensive   wet        34.854         9.091        11.950   reali     6     2 extensive   dry        50.387        19.245         0.000   reali     ...     26     0 intensive   dry        13.836        22.724         0.000 sinapis     27     0 extensive   wet         6.877         1.613        52.213 sinapis     28     0 intensive   wet        22.770         0.537        44.901 sinapis     29     0 intensive   dry        17.346        42.322         6.359 sinapis     30     1 extensive   wet        34.854         9.091        11.950 sinapis     31     1 extensive   dry        50.387        19.245         0.000 sinapis     ...      I'm interested in knowing if any of these variables influence the species and their respective counts. And this is the result of the "full" model.               glm(formula = count ~ type1 + type2 + perc.for.100m + perc.dry.100m +          perc.wet.100m + species, family = poisson, data = dej)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.8458  -1.1414  -0.4546   0.8297   2.2145            Coefficients:                    Estimate Std. Error z value Pr(>|z|)        (Intercept)    0.028129   0.523509   0.054  0.95715        type1intensive 0.196699   0.191960   1.025  0.30551        type2wet       0.071841   0.334286   0.215  0.82984        perc.for.100m  0.003741   0.008277   0.452  0.65130        perc.dry.100m  0.010952   0.010750   1.019  0.30829        perc.wet.100m  0.007467   0.011596   0.644  0.51960        speciessinapis 0.597837   0.187689   3.185  0.00145 **      Does this sound like the correct approach at all? _Some additional information_ As a side note, based on my exploration of the data, I would expect that the count would depend (at least) on the type2 variable, alas that's not what I got. ![enter image description here](http://i.stack.imgur.com/9sYc5.jpg) Using the "reverse logic", I tried if species can be predicted using my data, which presumably confirms the above results.               Call:     glm(formula = species ~ type1 + type2 + perc.for.100m + perc.dry.100m +          perc.wet.100m + count, family = binomial, data = dej)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.6322  -1.0136  -0.1568   1.0592   1.6407            Coefficients:                     Estimate Std. Error z value Pr(>|z|)       (Intercept)    -0.351192   1.658052  -0.212   0.8323       type1intensive -0.170583   0.651611  -0.262   0.7935       type2wet       -0.107377   1.078726  -0.100   0.9207       perc.for.100m  -0.002806   0.026807  -0.105   0.9166       perc.dry.100m  -0.010227   0.036982  -0.277   0.7821       perc.wet.100m  -0.006486   0.038071  -0.170   0.8647       count           0.345036   0.153811   2.243   0.0249 *     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 69.315  on 49  degrees of freedom     Residual deviance: 63.431  on 43  degrees of freedom     AIC: 77.431      _EDIT 1_ Aniko noticed that there may be an interaction between the type2 and species. Indeed!               Call:     glm(formula = count ~ type1 + type2 * species + perc.for.100m +          perc.dry.100m + perc.wet.100m, family = poisson, data = dej)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.0859  -1.1350  -0.1947   0.7109   2.7470            Coefficients:                              Estimate Std. Error z value Pr(>|z|)         (Intercept)             -0.357165   0.559987  -0.638  0.52360         type1intensive           0.196699   0.191960   1.025  0.30551         type2wet                 0.704769   0.429087   1.642  0.10049         speciessinapis           1.145132   0.306847   3.732  0.00019 ***     perc.for.100m            0.003741   0.008277   0.452  0.65130         perc.dry.100m            0.010952   0.010750   1.019  0.30829         perc.wet.100m            0.007467   0.011596   0.644  0.51960         type2wet:speciessinapis -0.962811   0.394038  -2.443  0.01455 *        _EDIT 2_ After removing the non significant terms (assuming I found the global maximum for the sake of data dredging) the story gets another twist into the right direction.               Call:     glm(formula = count ~ type2 * species, family = poisson, data = dej)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.7080  -1.1617  -0.1582   0.6979   3.1599            Coefficients:                             Estimate Std. Error z value Pr(>|z|)         (Intercept)               0.1542     0.2673   0.577  0.56408         type2wet                  0.6821     0.3237   2.107  0.03508 *       speciessinapis            1.1451     0.3068   3.732  0.00019 ***     type2wet:speciessinapis  -0.9628     0.3940  -2.443  0.01455 *


title: Classification in a applied financial setting

text: Many academic papers focusing on statistical learning in the applied setting of finance train a model such that their parameter set, B, defines the relationship y(t) = B*x(t) + e(t) y could be a 0-1 coded response for an increase (or not) of the asset price at time t. Ignoring the details of the model my question is as follows: given that you are predicting the response at the SAME time that you record your features, how would you make use of this model? Intuitively I think ok, my input space says we will see an increase in y but since they are at the same time, I can actually go and see what y has in fact done. Do I take the view of my model to be the truth and create a trading signal based on this information (for example)? Any practical insight would be greatly appreciated.


title: Am I using the correct model, if so what do I need to fix?

text: I am seeking to model the number of high tunnels on farms with 13 other factors. Factors include: lattitude, percent of the population with food access, percent of vegetable acres per county, proximity to the coast, percent direct sales, %of farms with a CSA program, %direct sales, etc. etc. I have been told to use the Poisson distribution and control for residuals, but the results are not great. Am I using the correct model? This model is intended to describe the conditions around high tunnels not to nessecarily predict where they are located. Currently I have an overdispersion factor of about 6.2...can that be improved?


title: How to convert a table into a graph in R

text: I have a table I want to convert into a graph (bar-graph or line-graph) The first column has fixed values. Twenty different values are simulated for these fixed values and kept in the next columns. I want to plot a graph of the fixed column against all the different twenty simulated columns. How do I go about it? I am using R2.12.1 > table name:bygrace V1 V2 V3 V4 V5   100 16 11 -6 1   120 -17 -12 7 -2   140 18 13 -8 3   150 -19 -14 9 -4   210 20 15 -10 -5                values for V1 are:100,120,140                   V2:16,-17,18,-19,20                 V3:11,-12,13,-14,15                 V4:-6,7,-8,9,10                 V5.1,-2,3,-4,-5      Actually, my table looks like the one above. The first column V1 represents Premiums charged by an insurance company,say last year (for 5 policyholders/policies) After critical analysis of the portfolio, the company decides to increase or decrease the premium amount for some of the customers. However, the company does not know exactly by how much it should increase of decrease the premiums for fear that it would lose customers. For this reason, four different scenarios/simulations are considered which are represented by the next four columns(V2 through V5) respectively. Now the task is to plot the premium amounts in V1 against these four different scenarios(bar-graph or line-graph; I think bar will be better). And my question is can this be done on one graph/at once? If yes, how shoul I go about it? Or do I have to plot the premium against each column separately? In fact, I have spent the whole of yesterday and today on this but I have not been able to get the desired result. Someone gave me an answer to try.And I am going to do that because it has given me an idea and I am very thankful to the one! Many thanks to everyone for their help Owusu Isaac


title: NaN p-value when using R's goodfit on binomial data

text: I am attempting to test the goodness of fit for a vector of count data to a binomial. To do so I am using the `goodfit()` function in the `vcd` package. When I run the function, however, it returns `NaN` for the p-value of the Chi- squared test. In my setup, I have a vector of count data with 75 elements.               > library(vcd)     > counts <- c(32, 35, 44, 35, 41, 33, 42, 49, 36, 41, 42, 45, 38, 43, 36,      35, 40, 40, 43, 34, 39, 31, 40, 39, 36, 37, 37, 37, 32, 48, 41,      32, 37, 36, 49, 37, 41, 36, 34, 37, 41, 32, 36, 36, 30, 33, 33,      42, 39, 36, 36, 29, 31, 41, 36, 39, 40, 37, 39, 39, 31, 39, 37,      40, 33, 41, 34, 46, 35, 41, 44, 38, 44, 34, 42)     > test.gof <- goodfit(counts, type="binomial",      +                     par=list(size=length(counts), prob=0.5))      Everything works fine, but when I inspect the `goodfit()` object I get the following:               > summary(test.gof)           Goodness-of-fit test for binomial distribution                                X^2 df  P(> X^2)     Pearson               NaN 75       NaN     Likelihood Ratio 21.48322 19 0.3107244     Warning message:     In summary.goodfit(test.gof) : Chi-squared approximation may be incorrect      I suspected it was a small sample size issue at first, but I also have a data set with 50 observations that does not return `NaN` for the p-value. I have also tried to switch the method in `goodfit()` to ML with similar results. Why would this function be producing `NaN` in this case? Is there an alternative function to calculate GOF on count data?


title: Is there a repeated measures aware version of the logrank test?

text: I am looking for a repeated measures version of the Logrank test. If I am correct, I am looking for an equivalent of the Friedman test for survival data. Any suggestions on where to look? (and R code will always be welcomed :) ) Thanks.


title: repeated measure ANOVA with non-consecutive repeats

text: I have a somewhat of a repeated measure setup but wanted to check if repeated measure 1-way anova can be applied here (assuming underlying assumptions of normality, sphericity etc are met). Situation: time series data of N subjects self-reporting a binary variable (yes/no), and also providing a continuous measure X with every self report. Hypothesis: measure X is different on days with "yes" reports compared to days with "No" reports. analysis idea: For each subject, separately average measure X for "yes" days and "no" days. X is a repeated measure for each subject. Run a 1-way repeated measure ANOVA. My main question is: if the repeated measure is not consecutively repeated (yes and no days can be in any sequence), seems like it shouldn't matter?


title: How to interpret lower order interaction when higher order interaction is significant?

text: I have a question about the interpretation of lower order interaction terms in the presence of a significant higher order interaction effect. Suppose I have a 2 (factor $A$) $\times$ 2 (factor $B$) $\times$ 2 (factor $C$) design where the highest order interaction ($A\times B\times C$) is significant and a lower order interaction term ($A\times B$) is also significant. Does the significant $A\times B\times C$ interaction render the $A\times B$ interaction uninterpretable (much like how main effects are rendered uninterpretable in the presence of a significant interaction)? Under this kind of circumstance, should I run a set of post-hoc/planned comparisons to check how different conditions are different?


title: Algorithm to Construct Covariance matrices

text: I am involved in a project where i need to construct some examples using non- negative Multivariate Gaussian Random variables(log-normal random variables). Part of the computation requires calculating the covariance matrices. Since i do not have any data and i am plucking numbers out of the air i was hoping if some one would be kind to verify that the following algorithm is valid or not.   1. I assume a bunch of positive log values as standard deviations for my variables. Lets call these $$s =[\sigma_1,\cdots, \sigma_n]$$   2. I then compute $$s^2 = s^{T}s.$$   3. I assume a correlation matrix $\rho$, which are all non-negative and this matrix is symmetric.    4. Multiply $\rho$ with $s^2$ element wise. This is a non standard operation, what i am suggesting here is the $(i, j)$th element of $s^2$ is multiplied by $(i, j)$th element of $\rho$. The matrix produced when all the element-wise operations are complete is the covariance matrix $\Sigma$. I have verified in mat lab that $\Sigma$ values of my examples are Choleski decomposable, as in i get a valid upper triangle matrix as result. Hence i am inferring from that they are positive-definite. These matrices are also symmetric. Can we infer based on the process and Choleski factorization, that covariance matrices are valid ? http://math.stackexchange.com/q/250912/23874


title: HMM a posteriori probs for hidden states and more

text: I have a two-state HMM and I am using Baum-Welch to estimate all the model parameters, including the Transition matrix. Then I use Viterbi to infer the optimal hidden state sequence. I am using such HMM to infer the hidden states, which are 1 or 0. However, currently the model gives me way more hidden states predicted as 1 than I'd like to have. What I'd like to have then is a good scoring system for the hidden states such that I can rank the states according to their probability of being 1/0 given all the data and parameters. I am thinking of two different approaches to limit the number of states predicted as 1 for my hidden states: **Q1.** Is it possible, and at which stage of my procedure, I should constrain the maximum number of states predicted as 1? Should I assume a certain form of the transition matrix, instead of estimating it, and then fit the other model parameters? I am open to suggestions here. **Q2.** I could also use some posterior probability of each state (given all data) to filter some states predicted as 1 out. How should I compute the posterior probability of each state predicted by Viterbi? I.e., I'd like to compute a vector of probabilities of the same length as the viterbi predicted sequence giving the a posteriori probability of each hidden state. Thanks.


title: Correcting for multiple testing on non-independent sliding windows

text: ## Question Values taken from adjecent windows in a sliding windows are corrolated. If I calculate a p-value from each window, how can I correct for the fact I have tested many windows, given that the windows aren't independent? ## Details I have data points (sequencing read start positions) that can take any one of a limited set of descrete values (genomic positions within a gene). Number of data points taking a given value is clearly binomially distributed: $$ X_i \sim B(n,1/L) $$ where n is the number of data points in total (the number of reads) and L is the number of values each one can take (i.e. the length of the gene). Thus the probability of getting $x$ or more reads starting at the same point is: $$ P(X_i \geq x) = 1 - P(X_i \leq x-1) $$ and the probability that any position in the gene has that many reads or greater is $ P(X \geq x) \times L $. Now consider the situation where we are intersted in the number of data points falling within a window of positions. The number of reads in a given window is clearly distributed thus: $$ W_i \sim B(n,s/L) $$ where s is the size of the window. We can easily calculate the probability that any given window will have $x$ or more reads. But what is the probability of _any_ window having $x$ or more reads. We can't simply multiply by the number of windows ($L-s)$ because the windows aren't independent, and this would rapidly lead to p-values of greater than one. So my question is, how do I calculate the probability that any window of size $s$ will contain $x$ or more reads?


title: Manually computing bootstrap standard errors in the linear regression setting

text: I have written an R script for obtaining bootstrapped standard errors in the linear regression setting. In practice, first in a model building step I select the final model to be applied at each bootstrapped sample (for simplicity suppose that it is a simple univariate linear model). Then I simulate `B = 1000` bootstrap samples (with replacement from the initial true dataset). For each simulated sample I fit the initially defined univariate regression model to the simulated data and store the estimated coefficient of the covariate of interest in `b_x`. Suppose now for simplicity that `B = 10`. At the end of the `B` simulations, I obtain the output matrix `boot_out` in this form:               boot_out <- read.table(text = "          iter      b_x         1       1.19         2       0.81         3       1.21         4       1.05         5       0.99         6       1.11         7       1.09         8       0.88         9       0.91         10      1.12",          header=TRUE)      Now I need to compute the bootstrapped standard error of the effect of the regressor of interest, call it `se_boot`. To this aim, I used:               se_boot <- sd(boot_out[,"b_x"])     se_boot      but I get a unexpected result: the bootstrapped standard error of `b_x` is **higher** than the standard error estimated through the initial model fitted to the true data.. By looking at my `b_x` values I found that the estimates of the partial effect of interest vary within a wide range, but I guessed that 1000 bootstrap replications were sufficient for refining the estimated standard error of `b_x`.. Before trying to do even more replications, could you please tell me whether the above passages are right? Thanks a lot for any help.


title: Multiple regression interaction with categorical IV with 3 levels

text: I need to carry out a hierarchical multiple regression. My data consists of one continuous dependent variable, 2 continuous predictor variables and a categorical IV with 3 levels.I need to assess whether an interaction effect occurs between each of the continuous predictor variables and the categorical IV in predicting my DV. I have read that it is necessary to dummy code k-1 variables with k being the number of levels in the categorical variable, is this correct in my case? If so then I would have 2 binary dummy variables. To create the interaction term do I then multiply my predictor variable x dummyvariable1 x dummyvariable2 and enter this term in the model? Or do I create 2 interaction terms: 1 created by multiplying my predictor x dummyvariable 1 and another by multiplying my predictor x dummyvariable2 and then enter both of these terms in the model simultaneously? I am an undergrad psychology student so I don't have a great knowledge of stats and have struggled to find an answer to this question in my reading that I can understand. I will be using SPSS. Any help would be really appreciated, even if you could direct me to some reading. Thanks!


title: Early split decision criteria for fast random (regression) forest estimation

text: Suppose I am on a node in a $regression$ tree and I am using running estimates of $\sum_{i \in Region_1} (y_i - mean(y_i)_{Region1})^2$ (and the same for Region 2) to determine whether to split the node at feature $X_j$ for large cardinality $|\\{X_j\\}| := p$, where $p/3$ is the number of randomly sampled features. My objective is to find some criteria to terminate the split decision early, perhaps at $p/5$ (for example), if the algorithm has stumbled onto a feature that gives a separation that is statistically very strong and unlikely to be non-trivially improved by further iterations of the splitting algorithm through the remainder of the sampled feature set. Have such approaches been tried? Are there any practical guidelines? How would I approach this in a statistically robust way (a way that does not wound generalisation - it should be a net benefit)? A quick way I could do it would be to generate a few trees on my data and plot the distribution of the split criteria calculation as we move from $1$ to $p/3$, for a given level in the tree, then use this to get my estimate of the probability that the feature at $p/5$ will not be beaten. But I am not sure how this would impact generalisation.


title: Reporting coefficient of determination using Spearman's rho

text: I have two non-normally distributed variables (positively-skewed, exhibiting ceiling effects). I would like to calculate the correlation coefficient between these two variables. Due to the non-normal distribution, I used Spearman's rank-order correlation, which returns a correlation coefficient and a significance (p) value. My results (n=400) show a significant ($p = 8 \times 10^{-5}$) but weak correlation (Spearman's $\rho$ = .20). If one uses Pearson's, one could describe the strength of the correlation in terms of shared variance (coefficient of determination, $R^2$ – in my case $R^2$ = .04, ie. 4%). Clearly, for Spearman's, it doesn't seem meaningful to square the $\rho$ value as Spearman ranks data. What is the best way to talk about the effect size using Spearman's $\rho$? Alternatively, following a discussion here (Pearson's or Spearman's correlation with non-normal data), I interpret the discussions as meaning that Pearson's correlation does not assume normality, but calculating p-values from the correlation coefficients does. Thus, I was wondering whether one could use Spearman to calculate the p-value of the correlation, and Pearson's to calculate the effect size, and thus continue to speak of shared variance between the two variables.


title: Find pattern in a signal - Variable Frequency and Amplitude

text: I assume this question is fairly simple but I haven't been able to find a way yet. I have a signal that I sample with 256 measures (might be 1024 as well) in live. The signal usually approaches a square like this: ![enter image description here](http://i.stack.imgur.com/XsX7w.png) My goal is to detect when the signal starts to not being square anymore and turns to a "Shark Fin" profile like this: ![enter image description here](http://i.stack.imgur.com/bDMFq.png) Ideally, I would have a function (filter?) that would compare the "Shark Fin" pattern ![enter image description here](http://i.stack.imgur.com/3y0v3.png) and the measured signal and output a value (in % for example) to express YES or NO the pattent is or isn't contained in the signal (with a certain % of similarity). The amplitude and the frequency of the signal are variable, which make it impossible to use a Matched Filter (please correct me if I'm wrong). I believe I need to look into Wavelet Transform but I couldn't understand how this works. Could you explain me how to apply this pattern detection in my signal? I haven't decided the language or the platform where the algorithm will be executed. The first step is to understand how the algo works. Thank you for your help.


title: Probability of x given past data and linear model assumption

text: I have a set of data where for every new point $x$ I need to know how probable this point was, given the past data. I can assume the data follows a linear model. I know I can find a simple linear regression and use $ \sigma^2 = \frac{1}{N-2}\sum{e_i^2} $ to get the variance of the error term. Then use that normal distribution to find the probability of my new $x$ value. Now here comes the difficult thing: But what if I can not assume that the model is correct and need to incorporate the uncertainity on slope and intercept?


title: mtry and unbalanced use of predictor variables in Random Forest

text: I am working on the Random Forest prediction, with the focus on the importance of predictor variables, and have a question regarding understanding of mtry and the actual usage of variables in the trees of Random Forest in R (package randomForest). I have a sample with over 1000 observations and a response vector with classification into 2 classes: 0 and 1. I have 4 independent variables. For the experiment I have chosen the following parameters:               rftry1=randomForest(x,y,xtest,ytest,mtry=4,ntree=500, importance=TRUE,keep.forest=TRUE,do.trace=TRUE,replace=TRUE,keep.inbag=TRUE,proximity=TRUE)      I then check which variables have been actually used in the forest, and find a very skewed usage:               > varUsed(rftry1)      [1] 2436 1758 1988 1156 The usage actually decreases dramatically from the first to the 4th variable. Would it actually affect prediction? Why does it happen like that? I have then rerun the experiment with mtry=2, and the usage has reduced overall but still is proportional to the usage of the first one.               MeanDecreaseAccuracy MeanDecreaseGini      1 0.006535855 2.177148 2 0.224706591 127.106268 3 0.006633846 5.020456 4 0.017522580 36.867821 I have noticed that the frequency of usage of the variable in the trees does not affect the importance to the extent of changing the order of importance of variables. So my main 2 questions are: 1) I guess I do not understand what exactly mtry controls ( I have read the number of variables sampled at each split). I have noticed that with mtry 2 and 4 some of the variables have been reused in the tree, while with m=1 they have not been and the trees were much shorter. Do i understand it correctly when I say: mtry=2 means that at every split instead of firmly choosing one next variable to split on we will randomly choose out of 2? If i have 4 variables overall, and variable 1 was used for very first split - what are my choices for left daugter and right daughter to split on? 2) Is mtry related to the unbalanced use of the variables for split? Why are the variables used like this, and does it affect the prediction outcome to make biased? All my variables are in the form of factors where 3 have 2 levels and one has 3 levels. Thank you for your answers!


I am confused with with the Poisson result and Exponential interpretation for spare part requirement analysis. I try to calculate the required number of spare parts for a disposable remove-replace item. From a real life case; if:   1. warranty period for a part is $t=400 hr$   2. **constant** failure rate for the part is $\lambda=0,003534$ 1/hour so mean time to failure $MTTF=1/\lambda=283{\rm hours}$ then if I use Poisson distribution to calculate what is the required number of spares NOT to be in shortage for 400 hours for this part with at least a probability of 95%, solution is 4 pieces of spares. (the reason is: 1st x value which satisfies below inequality is 4. $$ {P(x=0)\lor P(x=1)\lor P(x=2)\lor P(x=3)\lor P(x=4) }>0,95 $$ ) However logicly $MTTF$ says that average time between failure events will be 283 hours so even the part fails at 1st hour then I will need 1 piece of spare. So in this point I'm very confused since 1 and 4 are very different and has very serious consequences. (cost, customer satisfaction etc.) I also had the graph for Reliability for this part. (${\rm Reliability}=e^{-\lambda(t)}$ ) Observing the graph also didn't help me! Can you please explain "where and why my logic fails"? I think I've some trouble with understanding the true nature of exponential distribution and meaning of failure rate or MTTF. ![enter image description here](http://i.stack.imgur.com/jDNb7.jpg)


title: How to determine the quality of a multiclass classifier

text: Given   * a dataset with instances $x_i$ together with $N$ classes where every instance $x_i$ belongs exactly to one class $y_i$   * a multiclass classifier After the training and testing I basically have a table with the _true_ class $y_i$ and the _predicted_ class $a_i$ for every instance $x_i$ in the test set. So for every instance I have either a match ($y_i= a_i$) or a miss ($y_i\neq a_i$). How can I evaluate the quality of the match ? The issue is that some classes can have many members, i.e. many instances belong to it. Obviously if 50% of all data points belong to one class and my final classifier is 50% correct overall, I have gained nothing. I could have just as well made a trivial classifier which outputs that biggest class no matter what the input is. Is there a standard method to estimate the quality of a classifier based on the known testing set results of matches and hits for each class? Maybe it's even important to distinguish matching rates for each particular class? The simplest approach I can think of is to exclude the correct matches of the biggest class. What else?


title: Is expectation the same as mean?

text: I am doing ML at my university, and the professor mentioned the term Expectation (E), while he was trying to explain us some things on Gaussian processes. But from the way he explained it, I understood that E is the same as the mean μ. Did I understood right? If it is the same, then do you know why both symbols are used? Also I saw that E can be used as a function, like E($x^2$), but I didn't see that for μ. Can someone help me understand better the difference between the two?


title: t-distribution confidence intervals for non-Gaussian data but large n

text: I have a question concerning a claim I read in statistics books concerning the applicability of the t-distribution to compute confidence intervals for large $n$ if the data is not normally distributed (but has finite variance). The statement is that for small $n$ one should check the Gaussian assumption on the data, but if $n$ is large enough then the t-distribution can be used to compute the confidence intervals. The idea seems to be that the distribution of the mean will converge to a Normal distribution by the central limit theorem. Here is what I don't understand: A t-distributed RV can be characterized as a Normal RV divided by a $\chi$-distributed RV with $n-1$ degrees of freedom. I see that the distribution of the mean will converge to a Normal distribution, but what about the distribution of the standard error? If the original data is not Normal, why should the standard error follow a $\chi$-distribution? Therefore, why should the standardized RV have a t-distribution? The only way I could understand this is that for large $n$   * the distribution of the mean converges to a Normal distribution   * the t-distribution converges to a Normal distribution so basically for large $n$ one could use a Gaussian distribution anyway for the confidence interval. However, this somehow does not make the "detour" over the t-distribution. Or is there another mathematical reason why the mean divided by the standard error should have a t-distribution? Am I missing something? Thanks for your help.


title: Completely different results from lmer() and lme()

text: I got completely different results from lmer() and lme()! Just look at the coefficients' std.errors. Completely different in both cases. Why is that and which model is correct?               > mix1c = lmer(logInd ~ 0 + crit_i + Year:crit_i + (1 + Year|Taxon), data = datai)     > mix1d = lme(logInd ~ 0 + crit_i + Year:crit_i, random = ~ 1 + Year|Taxon, data = datai)     >      > summary(mix1d)     Linear mixed-effects model fit by REML      Data: datai             AIC      BIC    logLik       4727.606 4799.598 -2351.803          Random effects:      Formula: ~1 + Year | Taxon      Structure: General positive-definite, Log-Cholesky parametrization                 StdDev       Corr       (Intercept) 9.829727e-08 (Intr)     Year        3.248182e-04 0.619      Residual    4.933979e-01                 Fixed effects: logInd ~ 0 + crit_i + Year:crit_i                       Value Std.Error   DF   t-value p-value     crit_iA      29.053940  4.660176   99  6.234515  0.0000     crit_iF       0.184840  3.188341   99  0.057974  0.9539     crit_iU      12.340580  5.464541   99  2.258301  0.0261     crit_iW       5.324854  5.152019   99  1.033547  0.3039     crit_iA:Year -0.012272  0.002336 2881 -5.253846  0.0000     crit_iF:Year  0.002237  0.001598 2881  1.399542  0.1618     crit_iU:Year -0.003870  0.002739 2881 -1.412988  0.1578     crit_iW:Year -0.000305  0.002582 2881 -0.118278  0.9059      Correlation:                   crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y     crit_iF       0                                                   crit_iU       0      0                                            crit_iW       0      0      0                                     crit_iA:Year -1      0      0      0                              crit_iF:Year  0     -1      0      0      0                       crit_iU:Year  0      0     -1      0      0      0                crit_iW:Year  0      0      0     -1      0      0      0              Standardized Within-Group Residuals:             Min          Q1         Med          Q3         Max      -6.98370498 -0.39653580  0.02349353  0.43356564  5.15742550           Number of Observations: 2987     Number of Groups: 103      > summary(mix1c)     Linear mixed model fit by REML      Formula: logInd ~ 0 + crit_i + Year:crit_i + (1 + Year | Taxon)         Data: datai        AIC  BIC logLik deviance REMLdev      2961 3033  -1469     2893    2937     Random effects:      Groups   Name        Variance   Std.Dev. Corr         Taxon    (Intercept) 6.9112e+03 83.13360                       Year        1.7582e-03  0.04193 -1.000       Residual             1.2239e-01  0.34985             Number of obs: 2987, groups: Taxon, 103          Fixed effects:                    Estimate Std. Error t value     crit_iA      29.0539403 18.0295239   1.611     crit_iF       0.1848404 12.3352135   0.015     crit_iU      12.3405800 21.1414908   0.584     crit_iW       5.3248537 19.9323887   0.267     crit_iA:Year -0.0122717  0.0090916  -1.350     crit_iF:Year  0.0022365  0.0062202   0.360     crit_iU:Year -0.0038701  0.0106608  -0.363     crit_iW:Year -0.0003054  0.0100511  -0.030          Correlation of Fixed Effects:                 crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y     crit_iF      0.000                                               crit_iU      0.000  0.000                                        crit_iW      0.000  0.000  0.000                                 crit_iA:Yer -1.000  0.000  0.000  0.000                          crit_iF:Yer  0.000 -1.000  0.000  0.000  0.000                   crit_iU:Yer  0.000  0.000 -1.000  0.000  0.000  0.000            crit_iW:Yer  0.000  0.000  0.000 -1.000  0.000  0.000  0.000     >


title: Regression discontinuity design versus panel cointegration

text: I have a panel data (21 years) and I am trying to figure out whether I should use regression discontinuity design (RDD) or panel cointegration. I do have a randomness in the assignment variable and so I don't have any problem using the RDD; I can examine the causality (local average treatment effect). My questions are: first, why one should use RDD over the panel cointegration; second how can the causality from the RDD be linked with the Granger causality resulting from the panel cointegration (assuming that there is the cointegration and hence it is possible to estimate the vector error correction model). As far as I understand, when you have assignment variable and that variable is random, you always use RDD (although it is possible to use panel cointegration), but if you don't have an assignment variable and that variable is not random or if you don't have an assignment variable, then you proceed with the panel cointegration. I would really appreciate if you could let me know whether my understanding is correct. Please suggest the papers if possible. Note that my dependent variable is in growth form (e.g. growth of wealth)in RDD, but if I have to use panel cointegration, I will have dependent variable in level form for long run and in growth form for short run. Example: Let's suppose, I have a binary variable (treat) and the dependent variable is y. There are two different approaches in RDD: non parametric and parametric approach. In parametric approach, we use all observations and then the coefficient on the treat is called as causal effect (we can also use other control variables) . In cointegration approach, we use all observations and we can examine the long run relationship between treat and y (plus other variables). However, the concern here is that variable treat is binary while other variables (dependent and control variables) are continuous. I test all continuous variables and found to be integrated of order 1. But, I am not sure whether it is to test for the unit root of variable treat (binary variable). Does it make any sense?. This is important before we decide to use cointegration since it should all varaibles should be integrated of order 1. How do I incorporate the varaible treat in the cointegration? How can I then compare the results of cointegration and RDD?


title: Calculating bootstrap confidence intervals on second derivatives of a GAM object

text: I am trying to calculate confidence intervals on the second derivatives of a GAM curve. After estimating the derivatives, the idea is to then use function `boot()` to get confidence intervals, but I am unsure where to apply it along the way, or exactly what statistic(s) to bootstrap. The following code calculates second derivatives for specific values of `epi.i`               library(mgcv)     library(boot)          epi.i<-c(7, 8, 9, 10, 8.01, 6, 7.01, 5.999, 4.003, 12,      2.01, 7.02, 10.01, 8.2, 5.9, 3.9, 6.999, 4.0001, 3.99, 6.001, 8.001, 5.99,      7.9, 6.99, 3.98)          Rare25<-c(4.471429, 4.551474, 4.204894, 4.456710, 3.504348, 4.175824, 4.298193,     4.406838, 3.058707, 4.451128, 1.000000, 4.327893, 3.580541, 4.082432, 3.630869,     4.352130, 3.919075, 3.795205, 3.380952, 2.993347, 3.775886, 3.766723, 3.852396,     3.923977, 4.308840)          epismooth<-cbind(Rare25, epi.i)          G2<-gam(Rare25 ~ s(epi.i))     #new data for prediction     newDF<-with(epismooth, data.frame(epi.i=unique(epi.i)))          #prediction for smoothed estimates at each unique epi.i value     #with standard error     B<-predict(G2, newDF, type='response', se.fit=TRUE)          #finite difference approach to derivative following     #example from ?predict.gam     eps<-1e-7                    #finite difference interval     X0<-predict(G2, newDF, type='lpmatrix')               newDFeps_p<-newDF + eps     X1<-predict(G2, newDFeps_p, type='lpmatrix')               # finite difference approximation of first derivative     # the design matrix     Xp<- (X0 - X1) / eps          #second derivative     newDFeps_m <- newDF - eps          X_1<-predict(G2, newDFeps_m, type='lpmatrix')     # design matrix for second derivative     Xpp <- (X1 + X_1 - 2*X0)/eps^2          fd_d2 <- Xpp %*% coef(G2)     cbind(epi.i, fd_d2)      Any advice on how to proceed with confidence intervals from here would be greatly appreciated. The following code is what I have tried, though it's incorrect:               bootX0 <- function(formula, data,  indices) {     d <- data[indices,] # allows boot to select sample     X0<-predict(G2, newDF, type='link', se.fit=TRUE)     return(mean(X0))     }           resultsX0<-boot(data=epismooth, statistic=deriv2,     R=1000, formula=Rare25~s(epi.i))     resultsX0     plot(resultsX0)               bootX_1<- function(formula, data,  indices) {     d <- data[indices,] # allows boot to select sample     X_1<-predict(G2, newDFeps_m, type='link')     return(mean(X_1))     }           resultsX_1<-boot(data=epismooth, statistic=deriv2,     R=1000, formula=Rare25~s(epi.i))     resultsX_1     plot(resultsX_1)               bootX1<- function(formula, data,  indices) {     d <- data[indices,] # allows boot to select sample     X1<-predict(G2, newDFeps_p, type='link')     return(mean(X1))     }           resultsX1<-boot(data=epismooth, statistic=deriv2,     R=1000, formula=Rare25~s(epi.i))     resultsX1     plot(resultsX1)                    Xpp <- resultsX1$t + resultsX_1$t - 2*resultsX0$t/eps^2          fd_d2 <- Xpp %*% coef(G2)     hist(fd_d2)


title: What is the correct sequence of subjects and sources need to be mastered to be a Machine Learning master?

text: I just started to mt MS program on Machine Learning subject. I am curious about the correct flow of the subjects have to be learnt to be proficient on ML subject. In addition can you suggest some resources about? Thanks...


title: Outliers causing Heywood case in CFA using MLM estimator in lavaan

text: I am trying to perform a CFA on data (10 indicators: $n=300$) that is severely non-normal but continuous (counts of a clinical behavior over a period of weeks): many cases are at zero, a fair few between 1 and 30, and a handful at values sometimes much higher (even $> 1000$). I am using lavaan in R with the MLM estimator to accommodate for this distribution. I am primarily comparing a single common factor model with a model containing three correlated factors. If I use all cases, I get warnings indicating a non-positive definite matrix, and negative error variance (i.e. Heywood case). Because of the severity of the outliers, in line with Bollen's (1987) recommendations, I explored whether the extreme values were the cause of the negative error variance. I removed all cases 3 standard deviations above the mean ($n=16$) and this led to admissible solutions when I re-ran my models, and sensible parameter estimates. Am I interpreting these errors appropriately, and does this sound like a reasonable procedure to have taken in these circumstances? On a related note, is there a limit to how high the Satorra-Bentler scaling correction factor can be before one should be concerned? I see values of around 4 in some cases (unsurprisingly given the non-normality of the data), although I cannot locate any literature that comments on this issue.


title: Gaussian state space forecasting with regression effects

text: The R package dlm implements filtering and smoothing (`dlmFilter` and `dlmSmooth`) for models with regression effects, but forecasting is not (yet) available for these models:               mod <- dlmModSeas(4)+dlmModReg(cbind(rnorm(100),rnorm(100)))     fi <- dlmFilter(rnorm(100),mod)     f <- dlmForecast(fi,nAhead=12)     Error in dlmForecast(fi, nAhead = 12):        dlmForecast only works with constant models      How can I do this in R? Thanks for your help!


title: How to calculate sample size for simulation in order to assert some level of goodness in my results?

text: I am a stats newbie, so apologies in advance if I'm asking a braindead question. I have searched for answers to my question, but I find that many of the topics are either too specific, or quickly go beyond what I currently comprehend. I have some simulation work that includes large datasets which become infeasible to simulate exhaustively. For the smallest of my datasets, an exhaustive run presents the following distribution of results from a total of 9180900 tests. Result/Frequency:   * 0 7183804   * 1 1887089   * 2 105296   * 3 4571   * 4 140 What the numbers _mean_ does not matter; what matters is that the larger datasets I have can stretch into billions of tests, and become far too time consuming to run. I need to constrain the workload. I feel I ought to be able to sample from the full set of tests to derive a distribution for the sample, and infer (within some bounds) that the results of an exhaustive simulation would exhibit roughly the same distribution. There is no bias inherent in the tests which are run, so uniformly randomly choosing inputs ought to provide a valid sample. What I do not yet understand is how I should go about selecting my sample size. In particular, the distribution exhibits a strange tail, and I fear that sampling too small will lose the lower frequencies. (The 140 occurrences of '4' account for only 0.0015% of the population!) So, my question is, what is the best way of calculating a sample size with which I can assert some level of goodness in my results? Or, am I asking the wrong question?


title: Statistical validation of RandomForest models

text: I am currently working on a RandomForest based prediction method using protein sequence data. I have generated two models first model (NF) using standard set of features and the second model (HF) using hybrid features. I have done Mathews Correlation Coefficient (MCC) and Accuracy calculation and the following are my results: **Model 1 (NF): Training Accuracy - 62.85% Testing Accuracy - 56.38 MCC - 0.1673** **Model 2 (HF): Training Accuracy - 60.34 Testing Accuracy - 61.78 MCC - 0.1856** The testing data is an independent dataset (means not included in the training data). Since there is a trade-off in accuracy and MCC between the models am confused about the prediction power of the models. Could you please share your thoughts on which model I should consider for further analysis? Apart from Accuracy and MCC is there any other measure that I should consider for validation? Thanks in advance.


title: How to assess overfitting?

text: This is a follow-up of the question I posted earlier. I am assessing the two RF models which are generated using two different set of features > NF - Test_Accuracy > Training accuracy (500 features) > > HF - Test_Accuracy < Training accuracy (125 features) Testing and Training is done using independent data sets and the accuracy is derived from the average of a 5-cross validation. The difference between the models are in the number of features. I am afraid there could be a possible overfitting in one of the model (It is not clear to me which model could be overfitting, because I have used independent dataset and k-cross validation on the datasets). I would like to know what are the standard methods (tools/libraries) which can be used to assess overfitting.


title: Adjusting axis properties in a timeseries graph

text: Consider the following code:               require(zoo)     data <- read.csv(file="summary.csv",sep=",",head=TRUE)     data = zoo(data$compressed, as.Date(data$date))     data <- aggregate(data, identity, tail, 1)     days = seq(start(data), end(data), "day")     data2 = na.locf(merge(data, zoo(,days)))          par(bty = 'n')     plot(data2,xlab='',ylab='entropy (bytes)')      How can one:   1. Adjust both the horizontal and vertical limits o the drawn axis to match the start and end of data (as an example, while the y value may vary between 20 and 1525, the axis shows 0 and 1500).   2. Increase the horizontal resolution (at least add some minor ticks) when the timeseries spans for a couple of years.  Using:               plot(data2,xlab='',ylab='entropy (bytes)', xaxs = 'i', yaxs = 'i')      I get: ![alt text](http://i.stack.imgur.com/dLc7l.png) Which is not exactly what I had in mind. The vertical scale doesn't start with the exat minimum (around 25), nor does it end with the exact maximum. The horizontal scale starts and ends in the middle. Hints?


title: Inserting small arrows (triangles) in bottom axis pointing up with R

text: I'm using the following code to insert vertical lines in a timeseries plot:               abline(v=as.Date("2007-09-19"),col="red",lty=2)     abline(v=as.Date("2008-03-10"),col="red",lty=2)     abline(v=as.Date("2008-04-28"),col="red",lty=2)     abline(v=as.Date("2009-08-07"),col="red",lty=2)      But I would rather prefer a small triangle pointing up in the bottom axis, along with a label. Hints?


title: What's the relation between hierarchical models, neural networks, graphical models, bayesian networks?

text: They all seem to represent random variables by the nodes and (in)dependence via the (possibly directed) edges. I'm esp interested in a bayesian's point- of-view.


title: What are the assumptions for applying a Tobit regression model?

text: My (very basic) knowledge of the Tobit regression model isn't from a class, like I would prefer. Instead, I have picked up pieces of information here and there through several Internet searches. My best guess at the assumptions for truncated regression are that they are very similar to the ordinary least squares (OLS) assumptions. I have no idea if that is correct, though. Hence my question: _What are the assumptions I should check for when performing Tobit regression?_ _**Note:_** _The original form of this question referred to truncated regression, which was not the model I was using or asking about. I have corrected the question._


title: Eigenvalue decomposition in matlab

text: I'm implementing PCA using eigenvalue decomposition in matlab. I know matlab has PCA implemented, but it helps me understand all the technicalities when I write code. I've been following the guidance from here, but I'm getting different results in comparison to built-in function princomp. Could anybody look at it and point me in the right direction. Here's the code:               function [mu, Ev, Val ] = pca(data)          % mu - mean image     % Ev - matrix whose columns are the eigenvectors corresponding to the eigen     % values Val      % Val - eigenvalues          if nargin ~= 1      error ('usage: [mu,E,Values] = pca_q1(data)');     end          mu = mean(data)';          nimages = size(data,2);          for i = 1:nimages      data(:,i) = data(:,i)-mu(i);     end          L = data'*data;     [Ev, Vals]  = eig(L);         [Ev,Vals] = sort(Ev,Vals);          % computing eigenvector of the real covariance matrix     Ev = data * Ev;          Val = diag(Vals);     Vals = Vals / (nimages - 1);          % normalize Ev to unit length     proper = 0;     for i = 1:nimages      Ev(:,i) = Ev(:,1)/norm(Ev(:,i));      if Vals(i) < 0.00001       Ev(:,i) = zeros(size(Ev,1),1);      else       proper = proper+1;      end;     end;          Ev = Ev(:,1:nimages);


title: How to improve a Fraud Classification Model?

text: I built a classification model (Logistic Regression) in order to classify data in Fraud or Not Fraud. This data is related with online CNP (Card Not Present) transactions and after choosing some parameters that seemed to be related with fraud, I tested the model. For this I used a training set of 225000 examples, and a test set of 75000. I conducted two different tests, in the first one I had 7 parameters and managed to obtain a 96% accuracy in classification. The problem is that the number of Fraud cases is much lower than the Not Fraud ones, and so regarding the Fraud cases I got an accuracy of only 11% while on the Not Fraud 90 and something %. In the second test I included more parameters, making a total of 15. Same training and test set size, and although the overall classification dropped to 92%, regarding the cases of Fraud I got an improvement to 30%, and Not Fraud still around the 90s%. I would like to maintain the overall classification accuracy around 90%, and improve the accuracy on Fraud cases to something like 65~75%, but I can't find more parameters that seem to be relevant, to include in my model and I am stuck as other than that, no more ideas come to mind... Can someone please give me some hints or ideas on what to try next in order to try to achieve these goals? Also, I have another doubt. Because the values of the parameters that I am using, have a very wide range, I applied Feature Scaling and Mean Normalization over them. I have 300000 example samples (Training set - 225000 and Test set - 75000). My question is, for each of these sets should I calculate the respective Average of each column and the Max - Min, in order to convert the values to a smaller scale, or should I calculate the average and max-min, based on the whole sample (the 300000 examples)?


title: Are there any conditions on the data in Fuzzy c-mean clustering?

text: I want to cluster a graph using the FCM algorithm, I used the adjacency matrix of the graph as the data, and the "Euclidean" distance as metric. The problem is that the adjacency matrix is full of zero's (depends actually on the degree of the node) and my network is a kinda huge more than 2000 node. For the results I get many small size clusters and some big size clusters, that some time contain about 50% of the population and it is not right. I think my problem is from the data representation. is there any conditions that must be satisfied apart from the FCM algorithm conditions. Thank you in advance.


title: How to handle empty groups in a repeated measurements two way ANOVA

text: Suppose you have N individuals you would like to test. You test at the same time for two categorical conditions where condition one could be A,B or C and condition 2 could be W,X,Y,Z. You test every individual n for every combination of these two conditions and the result is a dependent variable:   \--> So basically you have 3*4 = 12 groups you would like to compare using a two way repeated measurements ANOVa Suppose now that you have a minimal criteria you're results have to meet: If they don't meet this criteria there not included in the group: so the number of individuals in a group can be less then N and not every individual can be found in every other group.   \--> still ANOVA can be performed if checked for sphericity. Now it comes: But what if one of the groups is empty. Not due to mistakes but simply due to the fact that no good results can be obtained in that group. You're fully aware off that. But at this point performing a two way ANOVA isn't possible anymore it seems. Any suggestions on how to avoid this or how to turn this around? Thanks in advance


title: Spatial autocorrelation -- GLM, autocovariate, MEM (Moran's eigenvector mapping)

text: I am currently working on two marine species distribution modelling and also on their overlap distribution. For this I use a binomial logistic regression model (GLM) with response being, respectively, presence/absence species1, presence/absence species2, presence/absence overlap. My question concerns my first results concerning the overlap response (1 is coded when I have a presence of the two species and 0 when at least one species is absent). Explanatory variables are different environmental parameters (temperature, depth, ect.). After selecting a final model I examinated spatial autocorrelation (SAC) and found in the error residuals using Moran's I index (6 first lags are significately correlated spatial). I decided to apply two different methods: autocovariate and the Moran eigenvector mapping. I am using R (package `spded`, `spacemakeR` and `packfor`), Dormmann (2007) "Methods to account fo spatial autocorrelation in the analysis of species distributional data - a review" appendix and Borcard et al. (2011) Numerical Ecology with R p. 275-277. My question is related to the results these two methods (see Figure enclosed). It seems that my GLM and GLM+ac models are better than the model GLM+MEM (MEM: Moran's eigenvector maps). Isn't that weird? Following literature, it appeared to me that MEM should be the best method of the three. Do you have any clue why that can happen? Can a theoritical reason be able to explain this or should I investigate in the way I apply the functions? I already applied the MEM function using two different neighbourhood matrices (one leaving "empty" neighbours, while other leaving no 0-regions). Maybe I should add that my data are on large scale (statistical rectangle of the North Sea) and that a factor to identical position but at different time (20 years covered) has been added to spatially separate them. Any help or clue would be welcome. ![Accuracy parameters](http://i.stack.imgur.com/0aHFC.jpg) ![ROCPlot](http://i.stack.imgur.com/xNKax.jpg)


title: car racing game win factors

text: I play a car racing video game on my tablet computer and I have collected data from 224 races. In the game, two cars are matched up head-to-head and race over a quarter-mile. I would like to determine the relative impact of certain factors on the likelihood of winning the race. A data file is available at http://csr.datamustflow.com/csr_racing_times.csv . This file contains the columns listed below. In the file, a 'y' means 'yes/present' and 'n' means 'no/not present'. Looking forward to any feedback you can give me on the best way to go about this analysis. Edit: A new file with the race pairings is now available at http://csr.datamustflow.com/csr_times_paired.csv Data fields: Make_Model: name of the car E: Engine upgrade that can be purchased with in-game cash; if present, should lower race time N: Nitrous oxide upgrade: if present, should lower race time T: Tire upgrade; if present, should lower race time B: "Blogger" if you lose the race, you don't lose any "game points"; not expected to have any effect on race time win: whether the car won the race or not time: The race time over a quarter-mile; lower time = equals faster race me: whether it was me racing the car in that race or not PP: An integer that describes the overall performance of the car; a higher number means a better, faster car


title: Dice probability for Yahtzee large straight

text: I had a discussion with a friend of mine about what is the best course of action in the following Yahtzee scenario: The first throw with 5 dice results in the following values:               1, 3, 4, 5, X, with X = 3, 4 or 5      Based on probability, what will be the best decision to make to end up with a large straight (1-2-3-4-5 or 2-3-4-5-6):   1. Try to throw a 2 with the X die in the remaining two throws   2. Take the 1 and the X dice to throw a 1 and a 2 or a 2 and a 6 in the remaining two throws. I claim the first option has the best probability to succeed while my friend claims it is the second option. I tried calculating the probabilities, but it has been a while since I followed my probability classes, so I am not so sure. For the first option I believe the probability is 1/6 + (5/6 * 1/6) = 11/36. But what about the second option? What is the correct way to calculate this probability?


title: About normal distribution, solve function issue

text: Please help me about this question: $$X\sim N(13,4) $$ $$P ( 20-d < X < 20+d ) = 0.998 $$ I have trouble finding $d$.


title: Dimension of weight vectors in SVMs

text: For a given set of features (say with dimension a) and for a given set of labels (say m labels), how to relate the given features with the weight vector of the SVM in general? Will it be equal to product of the two dimensions always, or is it required to find the dimensions of the both the feature space and the label space and then find the resultant dimension of the weight vector?


title: Estimating effect of latent variable in regression

text: I believe that independent variables $X_1,X_2$ affect the dependent variable $Y$ through a latent variable $Z$ such that $$ \begin{align} Y &= \beta_0 + \beta_1Z \\\ Z &= \operatorname{Logit}^{-1}(\beta_2X_1 + \beta_3X_2) \\\ \\\ Y &= \beta_0 + \beta_1\operatorname{Logit}^{-1}(\beta_2X_1 + \beta_3X_2) \end{align} $$ Is it possible to estimate $\beta_2$ and $\beta_3$, given $Y$?


title: Estimating of point in the density of bimodal distribution

text: I'm trying to estimate a specific point in a bimodal distributions, It's the point that I indicated in my picture with the A. ![enter image description here](http://i.stack.imgur.com/YQhXE.png) As example we can generate a mixture of guassian and then try to determinate the point A(x,f(x)). This problem is I think a problem of optimization and estimation. My question is, though, how would I go about that in detail? Maybe I'm just missing a good tutorial somehow.


title: Rescaling for desired standard deviation

text: Apologies for what is probably a very basic question. I have looked around both here and in the usual places and haven't had any luck. I have read that there are at least two methods for linearly transforming data so that you can give your distribution a certain desired standard deviation. What are they and are there cases where you'd want to use one method rather than another? Just for concreteness's sake, let's say you have test scores from 0-50, a mean of 35 and sd of 10, and you wanted to rescale so the sd is 15.


title: ANCOVA with repeating covariate

text: I wonder if someone can help me with an ANCOVA covariate question? Perhaps you may know the answer to my problem. Subjects were assigned to a treatment or control group, and I measured a (continuous) dependent variable. There is also a covariate, and it is continuous. But the covariate is not independent for each subject. **Example** Subjects are assigned to two different kinds of teaching method (treatment and control). In each condition, they do 15 tasks e.g. homework problems, and they get grades for each (continuous dependent variable). But the homework problems have varying levels of difficulty (continuous covariate) and we wish to control for that (also study it). **Is it ok to just treat problem-difficulty as a covariate?** Put another way, each record in the data is a single subject on a single problem. But the covariate is not indepedent for each subject on each problem. Instead, it is a function only of the problem; all subjects have the same value for problem-difficulty for a given problem. Is it ok to just treat problem-difficulty as a covariate? Thank you very much for any insights.


title: Does the standard error assume the mean is the centre of the distribution?

text: In statistics we can calculate a mean of a sample and the standard error of the mean. Let's say the mean of a sample is 2 and the standard error 1.5. If we repeatedly sampled from this population, it would be expected that the average distance of the sample mean from the population mean is 1.5. Up until now I think I've subconsciously thought that 2 is assumed to be the centre of the distribution, or in other words the population mean. But am I right in saying that 2 is just as likely to be in the tails of the distribution as it is the centre of the distribution?


title: Term used to describe a particular sample

text: Suppose that we know that captured data follow a particular distribution. We have 10 sample. 9 of them are close to the supposed distribution but the last sample seems to be very distant, or event "abberant". What is the term to use to qualify this sample ? Maybe "outlier sample" ? Thank you


title: How to measure Accuracy

text: I want to measure the accuracy of my GPS Receiver module. The real coordinates are obtained from the Google Maps, and the actual received coordinates are the ones that the GPS receiver received. I tried searching for statistical tools that can somehow support my findings. I tried using Chi-square test in Microsoft Excel, but I don't know how to interpret the results. How do I interpret the result? Or what statistical tool is really fit for this?


title: What is the distribution of the sum of independent normal variables?

text: Just need to check the answer for the following question: **Question**   Suppose $X$ and $Y$ are two independent standard normal variables: $X$ ~ $N(0,1)$   $Y$ ~ $N(0,1)$   What is the distribution of $X + Y$ ? **My Working**   $X+Y$ ~ N($\mu_1 + \mu_2$, $\sqrt{\sigma_1^2 + \sigma_2^2})$   $X+Y$ ~ N($0 + 0$,$\sqrt{1^2+1^2}$   $X+Y$ ~ N($0$, $\sqrt{2}$)   Does this look correct?


title: How to compare two groups when there are many outcome variables?

text: Suppose we have two patient groups, healthy and ill patients, and we measure the chemical concentrations of multiple chemicals. Now we want to determine which chemical concentrations are different for the two groups. **Can we just apply a separate student T-test for each of the chemicals?** It seems that for a large enough number of chemicals, you will always get a significant result for one of them.


title: Optimal data structure for 2-way repeated measures ANOVA

text: This must be a very basic question, but I'm trying to do a 2-way repeated measures ANOVA (in Statistica and/or SPSS), and was wondering:   1. in what way to (re)structure my data so it is most flexible for further analyses?   2. what way to go in Statistica or SPSS? (univariate GLM?) In my case:   * Subjects obtain a score (1-100) in describing (30!!) animals   * They do this 3 times, using only conceptual (C), visual (V), or auditory (A) information.   * The resulting data is currently formatted as follows, with one subject a case:                          dog        cat      ...     Sbj  C  V  A    C  V  A    ...     s1   64 78 34   61 63 39   ...     s2   78 89 31   68 77 45   ...     ...  .. .. ..   .. .. ..   ...      So factors are animal {dog,cat,...} and source {C,V,A}. I have tried a 4 variable format as well:               Sbj Animal Source Value     s1 dog C 64     ...      but:   1. This breaks the one-subject-per-case rule that the repeated measures ANOVA procedures in SPSS and Statistica rely on. In either case, SPSS requires me to specify all factor-levels manually which does not work as they are not variables, and moreover involves 20 levels! (in case of animal).


title: How do I run multiple linear regression with a limited data set for each subject?

text: I'm trying to use limited data across a range of variables to make predictions. There are ten variables and each subject has three of the ten variables defined. It's approximately random which three variables each subject has defined, but not completely. If each subject had all ten variables defined I'd know easily how to run a least squares linear regression. I have one subject who has all ten variables defined - using this data, I could try to fill in what the other seven variables would be for each subject, but it doesn't seem right. What's the best way to calculate the coefficients for a linear regression model in the case where all subjects only have three variables of the ten? Thanks!


title: Error propagation calculation yielding negative variance

text: I am trying to calculate the standard deviation of the sum X = A + B. A and B are mean values, and I do not have access to the source data. A is 0.46 with an SD of 0.014 (SDa) and B is 0.375 with an SD of 0.018 (SDb). Using SPSS and a different dataset measuring A and B in a different population I have concluded that A and B are correlated with a correlation coefficient of -0.4 and covariance of -0.00037 (COV). I find it propable that the relationship between A and B should be similar between populations and I am therefore then using this covariance as an estimation of the covariance in the dataset from where A and B are derived. To the best of my knowledge, the equation to calculate the SD of X = A + B (SDx) should then be: SDx = SQRT(SDa^2 + SDb^2 + 2*COV) This, however, gives a negative number to take the sqare root from, i.e. a negative variance for X, resulting in SDx being a complex number. A negative variance should not be possible, so I'm clearly making a mistake somewhere, but I can't figure it out. I couldn't find this specific question answered here before, so I thought I'd reach out and ask for some help. Any suggestions would be much appreciated! Best regards, Mattias


I have read about controversies regarding hypothesis testing with some commentators suggesting that hypothesis testing should not be used. Some commentators suggest that _confidence intervals_ should be used instead.   * What is the difference between confidence intervals and hypothesis testing? Explanation with reference and examples would be appreciated.


title: How to equalize histograms

text: I am learning some image processing stuff and equalizing a histogram comes up as an important topic. I have followed the procedure listed on Wikipedia but the resulting equalized histogram does not look much better to me. Let's say I have the following histogram where values only range from 0 to 7:               0 | 2     1 | 6     2 | 12     3 | 13     4 | 6     5 | 11     6 | 7     7 | 7      The equation to equalize the histogram is as follows: $h(v) = round(\frac{cdf(v) - cdf_{min}}{(M*N) - cdf_{min}} * (L - 1))$ So, using this I get the following: $h(0) = 0$ $h(1) = 1$ $h(2) = 2$ $h(3) = 4$ $h(4) = 4$ $h(5) = 5$ $h(6) = 6$ $h(7) = 7$ Transferring this over to the new histogram, I get the following equalized histogram:               0 | 2     1 | 6     2 | 12     3 | 0     4 | (13 + 6) = 19     5 | 11     6 | 7     7 | 7      This doesn't look right at all... Have I done something incorrectly?


title: Frequentism and priors

text: Robby McKilliam says in a comment to this post: > It should be pointed out that, from the frequentists point of view, there is > no reason that you can't incorporate the prior knowledge into the model. In > this sense, the frequentist view is simpler, you only have a model and some > data. There is no need to separate the prior information from the model Also, here, @jbowman says that frequentists use regularization by a cost/penalty function, while bayesians can make this a prior: > Frequentists realized that regularization was good, and use it quite > commonly these days - and Bayesian priors can be easily interpreted as > regularization. So, my question is, can frequentists in general incorporate into their models what Bayesians specify as priors? Taking the regularization as an example, is the cost/penalty function really integrated into the model, or is this a purely artificial means of adjusting the solution (as well as making it unique)?


title: Uncertainty of extrapolation (curve fitting)

text: I have estimated (MC simulated) some probability values **y** , that each depends on a value of **x** between 0 and 1. Say, for instance, that the vector **x** contains $x_{1} = 0.1,\ \ x_2 = 0.2,\ \ x_3 = 0.3,\ \ x_4=0.4,\ \ x_5 = 0.5,\ \ x_6 = 0.6,\ \ x_7 = 0.7$ and the vector of estimated mean values **y** contains $y_{1}(x_1) = 0.340,\ \ y_2(x_2) = 0.329,\ \ y_3(x_3) = 0.322,\ \ y_4(x_4)=0.299,\ \ y_5(x_5) = 0.278,\ \ y_6(x_6) = 0.255,\ \ y_7(x_7) = 0.237.$ I also have an estimated standard deviation value for each of the estimated values in **y**. Let these uncertainty values be contained in a vector **s**. Now I estimate $y(1)$ based on extrapolation (curve fitting). I do this by fitting my values to a curve on the form $y(x) = q \exp \left\lbrace -a\left(x-b\right)^c \right\rbrace$ by minimizing a weighted mean square error function (with weights according to **s** ) in order to find the optimal parameters $q, a, b, c$. Then I have obtained my estimate $y(1)$. However, my question is how to measure the uncertainty in this estimate. Is there for instance any good way to exploit the values in **s**? So far, I have thought of this way: The values in **s** was actually estimated by $s_i(x_i) = \sqrt{\frac{1-y_i(x_i)}{Ny_i(x_i)}}$ where $N$ was the sample size in the MC simulation of $y_i(x_i)$. Then I thought I could perhaps simply use $s(1) = \sqrt{\frac{1-y(1)}{Ny(1)}}$ as a measure of the uncertainty in $y(1)$, even though I never actually estimated $y(1)$ using $N$ samples in that way. Also, it's a drawback that this value doesn't depends on the size of **y**. Any suggestions? It would be great to obtain some sort of standard deviation value s(1), so that I can easily find the relative error $RE = s(1)/y(1)$ as well. (I have some other methods that I use to estimate the probability y(1), and it would be nice to compare $s(1)$ and/or $RE$ of the various methods...) Thank you! :-)


title: How to compare a data point with the k-means cluster

text: Say I have the centroids calculated via kMeans, and clusters identified. How do I take an incoming datapoint ( a row in matrix with all the features) to compare with the points and figure out which cluster the datapoint belong to? How to design a runtime like this? I assume we only have to calculate the clustering every x hours, and make sure of the cluster result every now and then. Any response will be appreciated! Machine Learning Newbie Simon


title: Consistent estimator for probability density function at a given point

text: I am trying to conduct a simulation study on a variation of the kernel density estimator. In that experiment, there is a parameter in one of my formulas that involves the value of the unknown pdf at a given point x. Therefore I need to find a consistent estimator to estimate the value of f(x), but I have no clues on where I should get started with. Can anybody suggest an estimator of f(x) that is consistent or which family of estimators I should be looking into, particularly, for a normal distribution? Thank you!!