title: What are Effective Regression Techniques for Linguistic Analysis of Linked Data?

text: Cross-post from MathOverflow where it was suggested that I might get better results here. I am in the early stages of a problem that involves parsing a large number ($\approx 5 \times 10^9$) of documents (web pages) and estimating values from them. In particular I need to identify pages that offer downloads of pirated files and estimate their traffic. Whereas some websites make their traffic statistics public, many do not. I am left with a sparse data set with which to fill in a large number of blanks. The good news is that flagged pages are verified by a human, who extracts any page view data they can find. This means I don't have to worry about writing software to recognize or collect that data when it's present. My hopes are to apply my techniques from this problem to many others if I'm successful. As you can see this is actually a two-part question. There is a classification stage, where I estimate whether a site is likely to yield a pirated file, and an estimation stage, where I use regression techniques to estimate how many views the page might have. My (naive) approach would be to extract n-gram count vectors from each document (n-word-gram, not n-character-gram), which gives me a mapping from a symbolic set into a (high dimensional) integer set. It's a lot of dimensions, but I've done it before with some success. I could use this n-gram "profile" to apply logistic regression or angle comparison for classification, then linear regression for estimation. There are some really sophisticated machine learning algorithms out there, but I'm not sure what sort of performance evaluation techniques I could use on them in the estimation phase. However, I'm sure there's a better way to do this. I can't imagine how a linguistic profile of a web page has much correlation at all with its traffic. I believe that site traffic has a lot more to do with the pages that link to it, and I could imagine using graph theory to model the propagation of page views. However, unless the site is flagged for viewing that data won't be collected. Only a small amount of pages will have their page views recorded. I could perhaps write a regular expression to match common phrases used when expressing "how many views" they got, but I'm not willing to commit time writing some sort of complex classification algorithm for that particular feature. I'm also sure there's a large number of other measurements I could extract from a web page that I'm not thinking of. Really it's a problem of converting web page content and link topology into meaningful variables that can be used in regression analysis. So my questions to you are these:   * What techniques have proven effective for classification of documents?   * Where would I even begin to learn about modeling and estimating internet traffic? I'm sorry if these are not rigorous enough questions, but at this stage I need more of a lay of the land than advice on a particular theorem. These subjects are not my primary area of expertise.


title: why is there a huge difference existed in coefficient of determination obtained from 10-fold cross validation?

text: I'm using **gradient boosting regression model** (GBRT). To evaluate this model, I use **10-fold cross validation** , in each of which I set same param and compute the coefficient of determination as a measure of fitting. However, I find that there exists a huge difference in coefficient of determination obtained from each fold, e.g., the coefficient of determination from `fold_1` to `fold_10` is:               [ 0.95310245  0.89725342  0.886711    **0.97063794**  0.84182142  0.80870443       0.70535911  0.8888032   **0.42510782**  0.70421155]      **Although the mean is 0.81 and std is 0.31, there is a fold in which the coefficient of determination is 0.4, while another fold is 0.97.** The only difference btw each fold is just the training & test data set, why does there exist such huge difference? Is such difference indicating that the performance of my model is not good?


title: Why is statistics useful when many things that matter are one shot things?

text: I don't know if it's just me, but I am very skeptical of statistics in general. I can understand it in dice games, poker games, etc. Very small, simple, mostly self-contained repeated games are fine. For example, a coin landing on its edge is small enough to accept the probability that landing heads or tails is ~50%. Playing a $10 game of poker aiming for a 95% win is fine. But what if your entire life savings + more is dependent on you hitting a win or not? How would knowing that you'd win in 95% of the time in that situation will help me at all? Expected value doesn't help much there. Other examples include a life-threatening surgery. How does that help knowing that it is 51% survival rate versus 99% survival rate given existing data? In both cases, I don't think it will matter to me what the doctor tells me, and I would go for it. If actual data is 75%, he might as well tell me (barring ethics and law), that there is a 99.99999% chance of survival so I'd feel better. In other words, existing data doesn't matter except binomially. Even then, it doesn't matter if there is a 99.99999% survival rate, if I end up dying from it. Also, earthquake probability. It doesn't matter if a strong earthquake happened every x (where x > 100) years on average. I have no idea if an earthquake will happen ever in my lifetime. So why is it even useful information? A less serious example, say, 100% of the places I've been to that I love are in the Americas, indifferent to 100% of the places I've been to in Europe, and hate 100% of the places that I have been to in Asia. Now, that by no means mean that I wouldn't find a place that I love in Asia on my next trip or hate in Europe or indifferent in America, just by the very nature that the statistics doesn't capture all of the information I need, and I probably can never capture all of the information I need, even if I have traveled to over x% of all of those continents. Just because there are unknowns in the 1-x% of those continents that I haven't been to. (Feel free to replace the 100% with any other percentage). I understand that there is no way to brute force everything and that you have to rely on statistics in many situations, but how can we believe that statistics are helpful in our one shot situation, especially when statistics basically do not extrapolate to outlier events? Any insights to get over my skepticism of statistics?


title: Why use Student's t distribution rather than Student's z distribution

text: So suppose I'm trying to decide whether to reject the null hypothesis that the mean of some random variable $X$ is zero based on the sample mean $\hat{\mu}$ computed from $n$ IID draws of this random variable. Then normally I would do a $t$ test $$t=\frac{\hat{\mu}}{\hat{\sigma}\;/ \sqrt{n}}$$ and check to see what the probability is that one would draw a more extreme value than $t$ from the $t$ distribution. But why even bother scaling by $\sqrt{n}$? Why not just use Student's z-distribution? After all isn't that the true distribution for the sample mean $\hat{\mu}=\frac 1n\sum_{i=1}^nX_i$? The only reason I can see to bother is so that for large $n$ Student's t distribution approximates the standard normal distribution, which might simplify computation, is this the only reason it's done?


title: Include Type B in effective degrees of freedom for expanded uncertainty?

text: This post has sat unanswered for a couple days in the TalkStats forum so I'm hoping somebody here can help me out. I'm looking for assistance in clarifying how Type B uncertainty is incorporated into the determination of expanded uncertainty. I'm using the Guide to Uncertainty in Measurement (GUM) as the primary reference. This lays out the steps for calculating expanded uncertainty in G.6.4. In particular, I need input on my interpretation of step 2, computation of the effective degrees of freedom (Veff) through the Welch-Satterthwaite formula (G.2b). From a plain reading of this step, I understand the combined uncertainty (Uc) to contain all Type A and all Type B contributions. In the denominator, the degrees of freedom for individual Type B components (Vi) is often set to infinity (as would be the case for a half-width from an equipment data sheet). This takes that term to zero, though the uncertainty component is still part of the combined uncertainty in the numerator. I'd be comfortable with this interpretation, except that upon reviewing an example on the NIST site I found that they were excluding the Type B components altogether when calculating Veff. Can anybody identify which technique is correct? Thanks for the help.


title: Find the data elements in a data frame that pass the rule for a node in a tree model?

text: So I have used the rpart package to create a tree model and I found an interesting rule and wondered if there was an easy way to see which observations in that data frame pass that rule. It seems very tedious to use path.rpart to find the path it took down the tree, and manually enter those filters into the data frame to look for them. Is there a method where I can pass a tree and/or a node, and a data frame and return all the elements in that frame that ended at that node?


title: Prediction from gam model at X = 0 not equal to intercept

text: I'm having a problem with GAM. Specifically when I try to make a prediction at a new data value at the origin, my prediction is not equal to the intercept. Can someone help me understand why this is happening? Here is an example:               set.seed(42)     N = 1000     x = runif(N)     y = x^.22+rnorm(N)     m = gam(y~s(x))     summary(m)     predict(m,newdata = data.frame(x=0))          Family: gaussian      Link function: identity           Formula:     y ~ s(x)          Parametric coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.78115    0.03164   24.69   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Approximate significance of smooth terms:            edf Ref.df     F  p-value         s(x) 2.158  2.687 13.79 7.61e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          R-sq.(adj) =  0.0344   Deviance explained = 3.65%     GCV score = 1.0043  Scale est. = 1.0011    n = 1000     1> predict(m,newdata = data.frame(x=0))             1      0.3756038       I think the prediction should be equal to the intercept when x = 0, but it's not!


title: Capture Recapture Method for Estimating Crowds

text: Does anyone know of a study that has actually used capture recapture methodology to estimate visitor numbers at an open air event? We are thinking of using this method at our annual festival, which has no fixed entry or exit point.


title: Difference between random effect and random intercept model

text: I am looking at clustered data and because I was trained in economics I tend to look at fixed effects and random effects as solutions. An alternative would clearly be multi-level modelling. However, for me the clustering is just a nuisance, so I only care about "controlling it out". However, looking at the equation of a random intercept model, I cannot see the difference to my random effect model. So I was wondering whether there just might be different names for the same thing, due to differences between disciplines? But then why can you (in Stata) estimate them differently?   * RIM: `xtreg y x, i(id) mle`   * RE: `xtreg y x, re` However, the equation for both seems to be: $$y_{ij}= b_{0} + b_{1}x1_{ij}+ u_{j} + e_{ij}$$ (Don't know how to format it nicely, but Google gives examples of RE and RIM equations, which look the same to me easily:   * RE: http://en.wikipedia.org/wiki/Random_effects_model   * RIM: http://www.bristol.ac.uk/cmm/learning/videos/random-intercepts.html


title: Iterative process for removing extreme samples

text: My samples follow heavy tail distributions. I use a process to detect and remove "extreme" samples that goes like this:   1. Measure mean and standard deviation of samples.    2. Remove samples higher than mean plus 4 standard deviations.   3. Repeat from Step 1 - for a total of 3 times. If there are no extreme samples then the above process removes none; if there are many then they skew mean/std.dev, however, the multiple iterations take care of that. For my problem, the above algorithm removes 0-5% of samples and empirical tests show that it works quite well and seems quite robust. However, is this process sound? Does it have a formal name to look it up? Notice that, I can not just find the "optimal" number of standard deviations above which we reject samples, since this whole process needs to be automated for multiple datasets and run as part of a live system. The datasets are similar but not identical: sometimes they are small and there are no outliers (no luxury to lose samples); sometimes they are large and the number of extreme examples that need to be removed varies and can be quite high. Data represent user actions (one variable per action type) and we have multiple sets of users we need to compare (different segmentations). However, many "users" are actually bots that typically (but not always) do the same actions multiple times (way more times than humans would do). It is fair to say that we have a mixture of two distributions that we can't really distinguish, human-users and bots. We care about differences of behaviours of human-users after we remove as many bots as possible. In the long term we will build a classifier for bots but for now we are after some quick solution.


title: Understanding Singular Value Decomposition in the context of LSI

text: My question is generally on Singular Value Decomposition (SVD), and particularly on Latent Semantic Indexing (LSI). Say, I have $ A_{word \times document} $ that contains frequencies of 5 words for 7 documents.               A =  matrix(data=c(2,0,8,6,0,3,1,                        1,6,0,1,7,0,1,                        5,0,7,4,0,5,6,                        7,0,8,5,0,8,5,                        0,10,0,0,7,0,0), ncol=7, byrow=TRUE)     rownames(A) <- c('doctor','car','nurse','hospital','wheel')      I get the matrix factorization for $A$ by using SVD: $A = U \cdot D \cdot V^T $.               s = svd(A)     D = diag(s$d) # singular value matrix         S = diag(s$d^0.5 ) # diag matrix with square roots of singular values.      In 1 and 2, it is stated that: $WordSim = U \cdot S$ gives the **word similarity matrix** , where the rows of $WordSim $ represent different words. `WordSim = s$u %*% S` $DocSim= D \cdot V^T$ gives the **document similarity matrix** where the columns of $DocSim$ represent different documents. `DocSim = S %*% t(s$v)` **Questions:**   1. Algebraically, why are $WordSim$ and $DocSimS$ word/document similarity matrices? Is there an intuitive explanation?   2. Based on the R example given, can we make any intuitive word count / similarity observations by just looking at $WordSim$ and $DocSim$ (without using cosine similarity or correlation coefficient between rows / columns)?  ![enter image description here](http://i.stack.imgur.com/3KJuT.png)


title: Multiple imputation with additive constraint

text: Say I have a dataset which includes three binary (dichotomous) variables: A, B, C, and C takes the value of A OR B. Many observations have missing A and B, and non-missing C. I guess that C should be included in any multiple imputation model, since it contains valuable information. But is there any multiple imputation package (preferably in R) giving the ability to specify the relation C = A OR B, so that only sensible values (satisfying the relation) are imputed? Thanks in advance!


title: Kalman Filter Correction efficiency

text: I was wondering if Kalman Filter used in a way to correct and reduce forecast errors is useful in real life forecast.Since we are using output forecast data and measurement data from t-1 to correct data at t. In real life we do not have the real measurements value so how will we be able to use Kalman Filter correction on our forecast values. Suggestion: What do you think of using at some point (t+n) Kalman correction of our forecast(at t+n) as the real measurement in order to use it to correct forecast data at (t+n+1)


title: cluster analysis, Ward: how to evaluate number of clusters and their quality?

text: I have a table of similarities (cosines) and I clustered it with the Ward method. Great outcomes, a wonderful dendogram, but then I tried to evaluate the quality of this cluster solution and I got stuck. First: identifying the number of clusters in my data (cause in Ward is not like k-means where you have to set a precise number of clusters). I calculated the sum of squares (see attachment) to see how many clusters are there, but there isn't a proper "elbow" in the data, so how many clusters shall I consider? Second: trying to calculate the purity of the clustering (with the tool CluTo), by indicating 4, 5, 6, 7... clusters, I can see that the purity increases the more clusters I indicate. Of course. If the number of clusters equals the number of instances of my data, then purity is 1 (the maximum). dah. Any suggestion on how to report this? (number of clusters? quality of the clustering solution?) ![enter image description here](http://i.stack.imgur.com/yI0Tq.jpg)


title: Determining best approximator based on repeated measurements

text: I have a set of data (which I have attached a plot of for your convenience) in which I have 2 categories of interest - let's call them A and B. For each of them I want to select from a series of categories (let's call them C_x) the one which best appproximates each one of them in terms of the measure at hand (in my case reaction times). ![enter image description here](http://i.stack.imgur.com/oHpgC.png) I have thought of 2 approaches:   * Take the one with the closest mean and perform a t-test. If it's significant - bingo!   * take my category of interest and run t-tests against all C_x categories and choose the one with the highest significance. I would like to know which of the above you think is best. Also, recently I have been told that my approaches have just a whiff of black-hat data mining. I fail to see why, but if this is the case could you suggest a better way of going about this?


title: Two independent events affect probability of third event: find CI for third event

text: I know the title is vague: I couldn't think of a good one. The basic question can be worded several different ways, but here is how I think of it to make it easier for myself: > Sally has only one umbrella. Every morning there is 50% chance of rain, and > every evening there is 50% chance of rain. Sally goes to and from school > every day, and if it's not raining she will not take her umbrella with her. > How many times does she have to walk in the rain without her umbrella? So the way I'm going about this is to generate an Excel model that tests at least 30 samples, with each sample having 30 back and forth trips, using some formula to determine whether or not that trip is made without an umbrella even if it's raining. I've noticed that whenever it rains (let this be 1), if the total previous number of 0s (trips where it doesn't rain) up to the previous 1 is odd, then Sally is screwed. For example: if we have 1 0 0 0 0 1, meaning that it rained on the first trip and the last trip, Sally is okay because the number of 0's (4) is even. If it's 1 then 1, the second 1 is 1 because I treat 0 as an even number. I'm having difficulty figuring out how to put this in some sort of Excel model. For example, I can have a column "# previous zeros," but how would I reset the formula once it hits a 1 in the other column? Also, how would I then proceed to take my results and probably construct a CI and five point summary in Excel? Thanks.


title: Bootstrapping estimates larger than full data estimate

text: In short, I am running bootstrapping on an association analysis data set in `R` (using the `Arules` package) and collecting the fit measures (support, confidence and lift) for the full data as well as the bootstrap samples. The number of bootstrap samples is equal 1000. The final resultant data set I have is the set of 1000 confidence measures for each rule as found by the Arules package. From this, I am calculating the 95% confidence intervals for the confidence measure of each rule. I noticed that for about 20% of the rules, the full data estimate of confidence is actually less than the lower 0.025th percentile of the set of bootstrapped confidence measures. This is a little counter-intuitive, since I would expect the full data estimate to be somewhere in the middle of the bootstrapped estimates (the bootstrapped sample being to the sample, as the sample is to the population etc.). What am I missing here?


title: Spectral analysis to find anomalies

text: With spectral analysis we are able to assess which series are smooth and which oscillate more quickly. Is it possible to find anomalies in a time series with spectral analysis? Could you please offer me a paper or reference about finding anomalies with spectral analysis?


title: Typical sequence of analyzing text

text: While going about examining textual data, do you have any process flow that works for you? What I do as of now:   * data collection   * importing into a corpus   * coding   * cleaning   * preprocessing: whitespace, lower case, punctuation, stopwords, stemming(maybe)   * frequency analysis   * associations, clustering or classification What should I be looking to do to improve from here?


title: How to check if two variables are dependent given a third variable?

text: In `R`, I have a dataset of many variables and based on correlation matrix I see that some of the variable are correlate with the others. For simplicity, let's assume that there are tree variables X, Y, Z that have high pairwise correlation coefficients. But it seems that information about the variable Y is already included in the variable Z. How can I check (in R or in Python) if X is dependent on Y given Z?


title: Does this autocorrelation plot look correct?

text: I have time series data that looks like this ![enter image description here](http://i.stack.imgur.com/0fc6a.png) and I'm getting an autocorrelation plot for this time series as ![enter image description here](http://i.stack.imgur.com/Ue1Ix.png) That is zero autocorrelation for every lag. But this would mean that the time- series is essentially random. Does this seem correct given the plot of the data? And if the data is essentially random, I suppose that it would not make any sense to construct an autoregressive model for it (i.e. if the data is random then data points should not have any dependence on its lagged values).


title: Sample size in BIC

text: The definition of the Bayesian information criterion is usually given as $BIC = -2 \text{ln}(L) + k\text{ln}(n)$, where $ln(L)$ is the maximized log- likelihood of the data given a particular model, $k$ is the number of free parameters that this model has, and finally $n$ is the number of data points that the model is being fit to. If I had a model which predicts the number of correct responses in two experimental conditions with 100 trials each, and my data consisted of the number of correct responses in each of the two conditions (e.g., 34 out of 100), then what exactly should I use for $n$, 2 or 200?


title: How to perform a meta regression with a random effect model?Which model should I use?How to start? (beginner)

text: I have to perform a meta-regression, using mixed or random effects model, but I don't have any software (except Matlab) and I'm new on this topic (having a relativelly poor statistics background). Briefly the problem I'm addressing is as follows. I've collected several effect sized from different studies. In my case, the effect size are percentages of adoption of a specific behaviour. For example, one study analysed a sample of 100 individuals and concluded that 35 % adopt the behaviour under analysis, whereas another study analysed a sample of 150 individuals and concluded that 45 % adopt the behaviour. My goal is to determine the average of adoption of the behaviour. All the studies are observational (no control groups). I have to perform a meta-analysis because its study have different characteristics, as for example, the sample size or for instance the average level of education of the sample. Some of the dependent variables are dummy variables. I read many things about meta analysis/regression but I'm so confused and I don't know how should I start my analysis, which model should I follow, what are the adequate formulation to calculate variance and the weights etc. How the heterogenity is included in the model? I thought about using this model: y=b0+b1x+b2x+...+u+e (fixed+random effect model) Is the way I'm thinking correct? Do you have any suggestion for me? If not, what is the step-by-step procedure? (in matlab or excel). I would appreciate if you could guide me on this analysis. I'm sorry if I wasn't clear or specific enough. Thank you in advange. Kind regards


title: interaction terms with opposite signs compared to the main variables

text: here I am again with another question. I am using interaction terms in my model. Problem is that no matter what specification or variable construct I try, interaction term always has opposite sign compare to my main variable. They both are significant as well. Any Idea why this is happening? is it due to possible multicollinearity? To add more detail, when I run the model without interactions i.e. $Y_t=α_1+β_1x_1+β_2x_2$ then I get right coefficients but when I add interaction of x1andx2 then in one model, main effect keeps its sign but interaction is negative and in another model main effect becomes negative and interaction becomes positive. Thank you.


title: Why and when create a R package?

text: I understand this question is quite a broad one, but I wonder what should be the decisive points in deciding to create (or not) a new package for R. To be more specific, I would add that the question is not about the reasons to use R in itself, more about the decision to compile various scripts and to integrate them in a new package. Amongst the points that could lead to these decisions, I have thought of (in a quite non-exhaustive fashion), of :   * the non-existence of other packages in the same sub-field ;   * the need for exchanging with other researchers and allowing reproducibility of experiments ; And amongst the points that could lead to a contrary decision :   * part of the methods used already present in some other packages;   * number of new functions not sufficient to justify to create a new independent package. I might have forgotten many points that could go in either list, and also, these criteria seem partly subjective. So, what would you say should justify, and at which point, to start bringing together various functions and data in a new documented and broadly available package ?


title: Variable for logistic regression is categorical and continuous so creates “missingness” in R

text: I am doing a logistic regression analysis using the glm command in R. It is to identify causes of valve narrowing beyond a certain threshold; 0=no narrowing, 1=narrowed. One of my variables is the size of a medical device that is implanted (range 25-36mm). Sometimes the device isn't implanted and I've left this as a blank field, but of course this is interpreted as a missing field. Not implanting the device seems to have a significant effect using Chi-sq analysis, and the size of the device has a significant effect using a t-test. How do I get around this in a linear regression model? To make it more complicated I actually have two different makes of the device: "C" and "D" with sizes 25-36mm, another device without a size "S" and then no device "N". Can it all be entered together or is it best to analyze separately outside of regression? What effect does the "missingness" have on various other variables that are in the analysis? Please & thankyou


title: Questions about the application of PCA

text: I'm learning about PCA. In general, I understand the concept and the underlying math, but I'm confused about a few things that I'm hoping someone can explain. Lets say that after performing a PCA analysis on some high dimensional data I find that the first three components capture 99% of the variance in my data. Now what? Is it possible to find which features comprise the first, second, and third components? It seems that a common approach is to plot the principle components. For example, see this page under `2.3.1. Dimensionality Reduction and visualization` for a PCA plot of the iris data. What is the interpretation of this plot? I see what I presume are two components: one component in red one component a mixture of green and blue? Is there a general "next step" to do after PCA. I understand the concept, but fail to see how it can be practically applied. Thanks for the help.


title: P-value NaN using chi2gof to validate lognormal distribution in dataset

text: I am trying to model a dataset of mine with a lognormal distribution using Matlab. I estimated the parameters via 'lognfit' and my generated datapoints with the fitted distribution look quite good compared to the observed data. To further verify, I wanted to use chi2gof to test the goodness-of-fit. Everything seems to work fine and the null hypothesis is not rejected but the p-value is NaN. I am using the approach applied to the the Weibull distribution under http://www.mathworks.de/de/help/stats/chi2gof.html So you can reproduce my result, this is an example that does not work.               >> test = lognrnd(-3.6299, 1.30985, 100, 1);     >> b = fitdist(test, 'Lognormal')          b =           LognormalDistribution          Lognormal distribution     mu = -3.62425   [-3.86888, -3.37962]     sigma =   1.2329   [1.08249, 1.43223]          >> [h, p, st] = chi2gof(test,'CDF', b)          h = 0          p = NaN          st =           chi2stat: 0.0183           df: 0        edges: [0.0015 0.2151 2.1372]            O: [94 6]            E: [93.6703 6.3297]      Also transforming the data to normal distributed datapoints and checking for a normal distribution results in a p-value of NaN. I am thankful for any help! Best regards EDIT: I was just able to make the transforming-to-normal approach to work: Would still be nice to know if this approach is valid and why the original approach does not work.               >> test = lognrnd(-3.6299, 1.30985, 100, 1);     >> test = log(test);     >> b = fitdist(test, 'normal')          b =           NormalDistribution          Normal distribution     mu = -3.68114   [-3.92861, -3.43368]     sigma =  1.24717   [1.09502, 1.44881]          >> [h, p, st] = chi2gof(test, 'CDF', b)          h = 0               p = 0.5330               st =           chi2stat: 4.1149     df: 5     edges: [-6.1930 -5.5713 -4.9497 -4.3281 -3.7064 -3.0848 -2.4632 -1.8415 0.0234]     O: [8 7 18 16 17 20 8 6]     E: [6.4814 8.9730 14.7440 18.9931 19.1820 15.1883 9.4282 7.0101]


title: Why did fixing the covariance value cause AMOS to mis-estimate variance?

text: I have 100 data points on two variables, `a` and `b`. The correlation between the two is .3 and the SD is 1. When I run the following in AMOS I get sensible estimates (.995 and 1.001) for the variances of both `a` and `b`, even though I've fixed the covariance to an inaccurate value. ![http://i.imgur.com/uYxZcal.png](http://i.stack.imgur.com/vHtj9.png) The estimates of variances stay about the same if I fix the covariance to be .3. However, if I fix the covariance to other values then the variance estimates can become wildly different. For instance, if I fix the covariance to .95 the variance estimates become 1.57 and 1.58 respectively. If I fix the covariance to .15 the variance estimates become .95 and .96 respectively. Why does this happen? In case it helps, here is the data in a Google spreadsheet.


title: How do principal components change upon the addition of new data?

text: How do the components from PCA change on addition of new data (i.e., $\frac{d(PC_1(x))}{d({\rm var}(x))}$? I am looking for any mathematical formulas and proofs (since it would be easy to understand).


title: Implement different variation of Holt Winter technique using R for 52 weeks of data

text: As per my job requirement I have to do forecasting using only Holt Winter technique in R.I have weekly data for 2 and half years & I have to predict weekly.I'm planning to build time series with frequency 52. I'm seeing that HoltWinters() function in R can do only Additive trend-Additive Seasonal (A,A) & Additive trend -Multiplicative Seasonal Model (A,M). I want to try other models like (M,A),(M,M). I saw these models can be done by ets() but the ets() has maximum frequency restriction of 24. Can you please suggest me some method in R where I can use these models with 52 weeks of frequency? Thanks in advance


title: Deriving the optimal value for the intercept term in SVM

text: I was reading andrew ng's machine learning lecture notes on SVM. I came across the following equation (finding the optimal value for the intercept term $b$ in the SVM problem): ![enter image description here](http://i.stack.imgur.com/0vD6t.png) However, I have no idea how the intercept term $b$ is derived by solving the primal problem ? I believe the primal's Lagrangian is: $$\min_{w,b} \max_{\alpha} \mathcal{L}(w,b,\alpha) = \min_{w,b} \max_{\alpha} \frac{1}{2} ||w||^2 - \sum_{i=1}^m \alpha_i [y_i (w^T x_i + b) - 1]$$ But how do I solve for $b$ ? Any help will be great. Thank you very much.


title: On which side resides the rejection area of a Chow test?

text: I know Chow test, it computes a test statistic, and the test statistic follows F distribution, but I don't know on which side does the null hypothesis reject.


title: Superiority test versus non-inferiority test

text: What is your best (or the best quote that you have) plain language explanation of the difference between tests for superiority and tests for non-inferiority? I _think_ a test of superiority a simply a garden variety two-sided hypothesis test for difference between two groups (e.g. treatments) aimed at providing evidence that a measure in one groups exceeds that in another; so H$_{0}\text{: }\theta = 0$, with _superiority_ being declared for a treatment if H$_{0}$ is rejected _and_ that treatment has greater effect. And I _think_ that a non-inferiority test is specifically a one-sided hypothesis test in which a new treatment is tested against an accepted standard treatment minus some amount ($\Delta$) which is accepted as being clinically irrelevant. So H$_{0}\text{: }\theta \le -\Delta$, with _non- inferiority_ being declared for the new treatment with respect to the standard if H$_{0}$ is rejected.


title: How can I include random effects into a randomForest

text: I'm not even sure that the question makes much sense, but I think I saw a couple of titles of papers where they proposed random forest with random effects. Is this possible in R?


title: Percentage outcomes in linear regressions

text: I have a study where many outcomes are represented like percentages and I'm using multiple linear regressions to asses the effect of some categorical variables on thes outcomes. I was wondering, since a linear regression assume that the outcome is a continuous distribution, are there methodological problems in applying such model to percentages, which are limited between 0 and 100?


title: Ploting classifiers' performance results for different datasets

text: I have five different data sets, and I applied five classifiers. I have also five performance evaluation metrics: precision, recall, F-measure, Accuracy, and AU-ROC. Is there any good plot that I can use to compare these results?


title: Data mining and time series : algorithm suggestion

text: I'd like to predict the behaviour (next action) of a internet user who had subscribed to a newsletter. 4 actions can be done by the user :   1. Don't open (/)   2. Read   3. Click   4. Buy  There is a hierarchy in this list. If you buy, you need to read and click. For the prediction, I have historical of the behaviour of the users : ![http://i.imgur.com/u25IP8U.png](http://i.stack.imgur.com/EbakD.png) As you can see, there is **3 difficulties** :   1. Two users can have two different subscription date (don't have the same origin in the timeline)   2. The users don't receive each email, because nobody subscribed to each subject (clothes, high-tech, sports)   3. I don't have any information about the subject of the email. For now I have in mind to use decision tree, logistic regression, svm,... For that I need :   * to transform the behaviour into a quantitative variable  Don't open (/) = 0 ; Read = 1 ; Click= 2 ; Buy= 3   * to pre-process the data and I have two main ideas : **1/ Create last behaviour variables like this :** ![http://i.stack.imgur.com/MaK8q.png](http://i.stack.imgur.com/D1esp.png) I loose a part of the time aspect here, because I don't know if there is a small or a big difference between the date of two actions. **2/ Create a score which include all the behaviour in one variable.** I can use a Gaussian kernel to had a decreasing weight with the time. ![i.imgur.com/JWaZNP3.png](http://i.stack.imgur.com/GDRmg.png) For this way other values than (0,1,2,3) can be tested ![i.stack.imgur.com/gdT5l.png](http://i.stack.imgur.com/hGDmu.png) Here I can't cross the behaviour of two consecutives behaviour in a decision tree for example. I didn't do that before, so if you have any suggestion, don't hesitate to purpose your idea.


title: How to perform discriminant analysis in R software?

text: I would like to perform discriminant analysis in R language. Please let me know the code and related packages for it.


title: Removing extraneous points near the centre of a QQ-plot

text: I'm trying to plot a QQ-plot with two data sets of about 1.2 million points, in R (using qqplot, and feeding the data into ggplot2). The calculation is easy enough, but the resulting graph is painfully slow to load, because there's so many points. I've tried linear approximation to reduce the number of points to 10000 (this is what the qqplot function does anyway, if one of your data sets is bigger than the other), but then you lose a lot of the detail in the tails. Most of the data points towards the centre are basically useless - they overlap so much that there's probably about 100 per pixel. Is there any simple way of removing data that is too close together, without loosing the more sparse data toward the tails?


title: Confirmatory factor analysis for data reduction (prior to regression)

text: I have a data set from an opinion survey with many variables, and to conduct regression analysis I would like to reduce the number of variables; because curently I actually have more variables than responses. I've tried principal component analysis(PCA)/exploratory factor analysis(EFA), but I find it difficult to "interpret" the different factors. However, I think one can group the variables well based on theoretical considerations alone, so I thought confirmatory factor analysis (CFA) might be an option. Unfortunately most resources on data reduction techniques point me to PCA/EFA. Does anybody know any resources that explain the process with CFA? Alterantively, can anybody help me with the following more specific questions? 1) If my sole purpose of is to use the factor scores in regression analysis, do i have to worry about the fit statistics of the CFA at all? It seems to me that all I need are the factor scores. 2) As far as I understand, one can either constrain the variance or set set one of the variables in each of the groups to 1, again for my purpose of data reduction, does it matter which approach I use? Does the choice affect my factor scores? Or more specifically does setting one factor loading in each factor to 1 the comparison across factor scores? 3) I am struggling with the actual computational step going from the CFA results to the factor scores. I was using the Lavaan package in R, and can't seem to get the factor scores. or do i have to calculate the manually using the loadings? Does anybody know how this is done, or can reccomend a different package? (I realise this is a more a programming question than a conceptual qustion, but thought I'll add it nonetheless, hoping that somebody who might know about CFA in theory also has practical experience with it in R) Thanks in advance! NB: In case anybody is alaremd by the statement about the relatively sample size above. I don't really care so much about inference here, as my sample is pretty much the entire population. NB2: I am planing to use the factors only on the right-had side of the regression. The left-hand side variable comes from elsewhere.


title: Algorithm for determining if strings are the same? (given the set)

text: Suppose I have a set of strings that are all related. (Imagine a list of ice cream flavors.) For any pair of strings, I want to determine the probability that they actually refer to the same thing. Is it possible to improve on simple algorithms (ie, Levenshtein distance) by taking the other members of the set into account?


title: Hoops in a Field

text: I have a simple question, but I haven't been able to find the answer I'm looking for anywhere. I'm sure there is some simple theorem or distribution that tackles this problem, but it evades me. Perhaps I haven't googled the right terms or looked in the right textbooks! My problem is as follows. Suppose a proportion ($p$) of a field is covered in water. The water is randomly dispersed in the field. The field has an area $A$. Now we fly over the field in a helicopter and randomly drop a hoop of size $a$ (where $aN = A$) into the field. The hoop lands at a random spot within the field borders. What is the probability that there will be water inside the hoop? My intuition about the problem is that for $a > A(1-p)$, the probability is 1 (because the area excluded by the hoop is smaller than the total area of the water), and as $a \rightarrow 0$, the probability approaches $p$. But what about for values of $a$ in the middle? Can the question be answered? I have a sneaking suspicion that the answer is embarrassingly simple. Thanks for your help!


title: Fitting t distribtution to financial data

text: The data I am using can be found here: http://uploadeasy.net/upload/nwv0.rar The variable is called "alvsloss". I want to fit a distribution to my financial data. First of all, I started with a normal distribution. The pdf is given by: $\frac{1}{\sqrt{2\pi\sigma^2}}\operatorname{exp}\left\\{-\frac{\left(x-\mu\right)^2}{2\sigma^2}\right\\}$ and the log likelihood, which is to be maximized is given by $ \ln\mathcal{L}(\mu,\sigma^2) = \sum_{i=1}^n \ln f(x_i;\,\mu,\sigma^2) = -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln\sigma^2 - \frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2. $ which gives as a solution $ \hat{\mu} = \overline{x} \equiv \frac{1}{n}\sum_{i=1}^n x_i, \qquad \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \overline{x})^2. $ I can calculate this in R via mean(alvsloss) and sd(alvsloss) and get the values -0.0004270872 and 0.02304159 or I do               library(MASS)     fitdistr(alvsloss, "normal")      and get the solution                      mean             sd             -0.0004270872    0.0230371244       ( 0.0004535429) ( 0.0003207033)      Ok, I understand this. But now my problem arises, I tried this with the t distribution. The pdf is given by: $\textstyle\frac{\Gamma \left(\frac{\nu+1}{2} \right)} {\sqrt{\nu\pi}\,\Gamma \left(\frac{\nu}{2} \right)} \left(1+\frac{x^2}{\nu} \right)^{-\frac{\nu+1}{2}}\\!$ the log likelihood can be found here: http://de.scribd.com/doc/52587265/80/STUDENT%E2%80%99S-t-DISTRIBUTION formula (17.8) (x is denoted by $\epsilon$. I am not understanding where the $h_i$ comes from in the formula, but this is not my main question.) First question: Seems to me, that there is no closed form solution? Second question: I do `fitdistr(alvsloss, "t")` and get the output                 m               s               df             -0.0004919768    0.0130128873    2.6340459185       ( 0.0003182568) ( 0.0003453702) ( 0.1620424078)     Es gab 28 Warnungen (Anzeige mit warnings())      I am not understanding, why I have three parameters now? The df is the $\nu$, but the others? What is R doing? Since there is no closed form solution, I guess it is using a iterated algorithm? Should I bother about the 28 warnings? They are: Warnmeldungen:               1: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     2: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     3: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     4: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     5: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     6: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     7: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     8: In log(s) : NaNs wurden erzeugt     9: In log(s) : NaNs wurden erzeugt     10: In log(s) : NaNs wurden erzeugt     11: In log(s) : NaNs wurden erzeugt     12: In log(s) : NaNs wurden erzeugt     13: In log(s) : NaNs wurden erzeugt     14: In log(s) : NaNs wurden erzeugt     15: In log(s) : NaNs wurden erzeugt     16: In log(s) : NaNs wurden erzeugt     17: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     18: In log(s) : NaNs wurden erzeugt     19: In log(s) : NaNs wurden erzeugt     20: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     21: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     22: In dt((x - m)/s, df, log = TRUE) : NaNs wurden erzeugt     23: In log(s) : NaNs wurden erzeugt     24: In log(s) : NaNs wurden erzeugt     25: In log(s) : NaNs wurden erzeugt     26: In log(s) : NaNs wurden erzeugt     27: In log(s) : NaNs wurden erzeugt     28: In log(s) : NaNs wurden erzeugt      And how can I plot this distribution with the estimated parameters?


title: Empirical Conditional Intensity Function

text: I'm trying to compute the empirical conditional intensity function in my data. I have two types of events, event Type1 and event Type2, that can occur in a time interval of, say, 1 day (events don't occur simultaneously). I'm assuming that when an even Type1 occurs the likeliness of another event Type1 occurring increases (self-excitation) and that when an event Type2 occurs the likeliness of an event Type1 occurring increases as well (cross-excitation). I want to figure out what happens to the empirical density of Type1 and Type2 events conditional to an event of Type1 occurring (in the specific, I want to see what happens to it during the 1st second after the occurrence of a Type1 event). If one looks at the non-conditional intensity than all you do is take my time interval (0,T), divide it into n non overlapping buckets and simply count the number of occurrences in each bucket. The way I'm doing it is to start a bucket of 1 second after the occurrence of each Type1 event, and divide the 1 second bucket into 20 sub-buckets of 0.05 seconds each. Then I count the number of Type1 and Type2 events occurring in every sub-bucket and repeat this for every bucket that I start after the occurrence of a Type1 event. Finally, I take the average number of occurrences for each of my sub-buckets across all the buckets and I plot my result. Is this the correct way of doing it? is it a problem that the buckets (the 1 second long ones that I start at every occurrence of a Type1 event, not the sub-buckets!!) overlap since I have Type1 events that occur more frequently than every 1 second? When I compare it to my model implied conditional intensity function the shape "fits" but the scale is totally wrong!


title: Joint pdf of a continuous and a discrete rv

text: Let us consider a manufacturing system. It involves 2 independent components. If one of these components fails then the entire system fails. Let $Y_j$ be distributed $\exp(Q_j)$ where $j=1, 2$. If component 1 fails first, then $Y_1$ is observed but $Y_2$ is not ($Y_2$ is censored). If component 2 fails first, then $Y_2$ is observed but Y_1 is not ($Y_1$ is censored). Therefore, if the system fails, we can only observe $u = \min(Y_1, Y_2)$ and the binary random variable $V$, which is $1$ if $Y_1 < Y_2$ and $0$ otherwise. How can I derive the joint pdf of a continuous variable $u = \min(Y_1, Y_2)$ and a discrete variable $V = 1$ if $Y_1 < Y_2$ and $0$ otherwise?


title: Applications of conditional covariance matrices

text: I have been trying to understand the methodology and application of conditional distributions in a multivariate normal setting as outlined in this Wikipedia page. In the context of that page, if we know the exact value of $x_2=a$, we can derive the conditional variance covariance matrix for $x_1$. Now, e.g., if I have 3 variables $x_1$,$x_2$ and $x_3$ in a MVN distribution, and I get the conditional distribution of $x_1$ and $x_2$ given that I know $x_3$ (either as a single value or as a distribution?), what information is this derived conditional distribution giving me? What are the applications of this distribution? Can I use this conditional distribution to make out-of- sample predictions or at least predictions at values within the range of the $x_3$ not present in the original data. I would appreciate if someone could point me to some examples of simple applications or literature. I'm not too technical but will do my best to read through and understand.


title: Normalising Constant for exponentiated function

text: What would the normalising constant be of the following, or atleast an approximation? I would like to avoid sampling. $$f(\theta)=\exp(-k_1e^{-k_2\theta^2}-\theta^2)\qquad\theta\in(-\infty,\infty), \quad k_1,k_2>0$$ In the event that its not possible, what would be a valid **polynomial** lower bound on $\log(f(\theta))$, that would take into account $k_1$ For anyone, who is interested, I end up with this expresssion because I place a Gaussian prior over (a squared exponential) kernel parameters, and trying to use variational bayes (hence, a lower bound).


title: Is it possible to compare multiple proportions in subdata sets?

text: There is one medical data set which contains 3 treatments: A, B, and C. The response variable is Yes/No. However, due to a specific medical reason, the dataset needs to be split into several sub-datasets. The major research question is to compare the difference of proportions between $P_A$, $P_B$ and $P_C$. If I did not need to split data, I could use a "simultaneous confidence interval" {Agresti et al., (2008)} to perform a multiple comparison between these three treatments. However, I can't find an appropriate statistical method to compare multiple proportions when a dataset is split into several sub-datasets. Can someone please give me a hint on this? Many thanks in advance.


title: Estimation of probability of a success in binomial distribution

text: Let's say we have two biased coins. The probability of tossing a head on the first coin is $\alpha$ and the probability of tossing a head on the second coin is $1-\alpha$. We toss both coins $n$ times and we say that there is a success when there is a head on both coins. If we denote this random variable by $X$ then $$X\sim B(n,\alpha-\alpha^2).$$ The question is how to properly estimate $\alpha$. The problem with a 'standard' approach ($\alpha - \alpha^2=\frac{k}{n}$ where $k$ is number of successes) is that the number of successes might be greater than $0.25n$ and in such a case we obtain that $\alpha$ is a complex number. We can treat such cases as $\alpha=0.5$, but then this estimator is not unbiased anymore.


title: Survival analysis for nested, censored and dependent data

text: I have to analyse the effects of different treatments on the survival of individuals for 1 week. But:   * the data are ' **grouped** ': individuals were grouped by 10, and the survival probability is affected by the number of individuals who died already   * the data are **nested** : I've got replicates.    * the data are **right-censored** : I want to include the number of individuals who didn't die after a week.  I use R, and I can't really find any packages to use for my problem. Cohort ad frailty analyses dont seem to ba appropriate. I was thinking maybe something from **epidemiology** for the first problem (and analyzing it as a 'contamination probability' or something).


title: How to create a multivariate Brownian Bridge

text: It is known, that a standard multivariate Brownian bridge $ y(\mathbf u) $ is a centered Gaussian process with covariance function $$ \mathbb E(y(\mathbf u) y(\mathbf v)) = \prod_{j=1}^d (u_j \wedge v_j) - \prod_{j=1}^d u_j v_j $$ I am not sure about how to constuct such a multivariate Brownian bridge. My first thought was to start somehow with a univariate Brownian bridge. I have found information about that and even a package in R that can do this, but only for the univariate Brownian bridge. I found this, but as I understand it, what has been done there is not a standard multivariate Brownian bridge as defined above or e.g. in this paper. I would appreciate any hints and support.


title: Understanding the conditioning in a GARCH process

text: In a GARCH model like the following $$y_t=\sigma_tz_t,\\\ \sigma_t^2=\omega(1-\alpha-\beta)+\alpha y_{t-1}^2+\beta \sigma_{t-1}^2$$ where $z_t$ is assumed to be iidN(0,1), we say that conditional on past information $y_t$ has the gaussian density $$f(y_t|y_{t-1},\sigma_{t-1}^2)=\frac{1}{\sqrt{2\pi\sigma^2_t}}exp\left(\frac{1}{2\sigma^2_t}y_t^2\right)$$ Am I correct in making the following conclusion?   * So when conditioning on the past, we know what the value of $\sigma^2_t$ is. Consequently, we can state that $y_t$ is conditionally normally distributed as $N(0,\omega(1-\alpha-\beta)+\alpha y_{t-1}^2+\beta \sigma_{t-1}^2)$.


title: Fisher's noncentral Multivariate hypergeometric distribution

text: I'm trying to calculate the probability of at least a given number of n ball colors in a urn where I have m > n balls. I would like to generate as many samples as I can of an urn with m balls, and then see how much of this balls are from the balls I'm interested (can be n balls < m). In R with `BiasedUrn package` I have density function and random number generator. Doing               rsamp <- rMFNCHypergeo(nran=1000000,              m=NumberBallsInUrn(vector > 50),              n=SampleNumber,              odds=rep(1, length(NumberBallsInUrn)))      I generate the sample I want, do you think that calculating the probability as               sum(rsamp[28,] >=12 & rsamp[26,] >=6)/1000000      could give me a fair probability of having at least (>=) 12 balls of ball 26 and at least 6 (>=) balls of ball 28 in the urn, similar to cumulative distribution of the hypergeometric distribution univariate case, where I calculate P(X > x) with               phyper(q, m, n, k, lower.tail=FALSE)      ? Could this be a fair calculation for balls drawn in an experiment in biology?


title: Can I use 3-way ANOVA for testing the difference between small proportions?

text: I'm looking at daily purchasing data over a period of several weeks. Specifically, we're interested in the % of users who make purchases in our app (it doesn't matter the type of purchase or amount spent for our purposes). On any given day, there might be ~10K active users and only about ~30 who actually make a purchase (it's a tough market!). There are three factors that we are interested, call them A, B, and C. Each factor has two levels. A few rows of data look like:               #   NU      NP  %buy    A   B   C   date     1   12407   24  0.0019  0   0   2   2013-03-06     2   11774   37  0.0031  1   0   2   2013-03-05     3   12331   39  0.0032  0   0   2   2013-03-04     4   12679   41  0.0032  1   0   1   2013-03-03     5   12555   61  0.0049  0   1   2   2013-03-02     6   12594   50  0.0040  1   1   2   2013-03-01     7   12466   71  0.0057  0   1   2   2013-02-28     8   12089   61  0.0050  0   0   1   2013-02-27      Here, each row is a daily observation of one group of users with a set of A,B,C factors, NU is the number of total users in that group, NP is the number of users in the group that made a purchase, and %buy is just the % of users in the group who made a purchase. For convenience, let's assume that there are no important time-varying trends in any of the groups, which obviously could complicate the analysis even more (and over a short time period that we're looking at is probably an ok assumption anyway). My first thought was to do 3-way ANOVA using the %buy as the observational outcome, but then I know I'm neglecting the importance of the sample sizes. Also, I'm thinking that the small proportions might mean a different test is necessary than run-of-the-mill ANOVA. **So my question is really what is the proper test that I should be running here to look for significant differences between groups and within groups?** Hopefully, I've given enough explanation to make this question clear.


title: How can I account for a nonlinear variable in a logistic regression?

text: How can I account for a nonlinear independent variable in a logistic regression? For example, consider this data set:               b1  b2  b3  b4  b5     1   1   0   0   20     0   1   0   0   20     1   0   0   0   13     1   1   1   0   8     1   1   0   1   5     0   1   0   1   4     1   1   0   1   5     0   0   0   0   8     0   0   1   0   13     1   0   0   1   20     1   0   0   0   29     1   1   0   1   40     1   0   0   0   53     0   1   1   0   68      Suppose I gathered some data like the above (suppose I had 10,000 rows) and I want to use them to predict future sets of similar data with an equation generated by a logistic regression. I get the above data. My independent variables are b1 through b5. By looking at the data, I can clearly see that independent variables b1 through b4 are binary (maybe they are male/female, rent/own, etc). However I see that b5 is clearly not binary. Upon graphing it and looking at the numbers I determine that the b5 variable has a bathtub shaped U curve with an equation of b5=(x-4)^2+4. I understand that using the b5 variable as it is would screw up my logistic regression - since any interval independent variables must have a linear relationship with the dependent variable. Yet, I must account for the variable b5 in my data set while doing a logistic regression for it to be accurate and predict well (right?). Given the above data, how would I go about making a predictive model with a logistic regression? Any references online or in books would greatly be appreciated as I am having trouble finding a good explanation of the above. Cheers =) EDIT: Sorry if my question is unclear. Consider this: Suppose we had the following data on many people: Y: Whether they had a stroke before the age of 50 or not 1\. Whether they own a pet or not. 2\. Whether they graduated college or not. 3\. Whether they are male or female. 4\. The number of hours they exercise every week. We want to train a logistic regression model based on historical data to predict if someone in the future with the same data will be more apt or not to have a stroke before the age of 50. The problem I am having is with the last variable. A continuous ratio variable like "number of hours exercised each week." Exercising more each week may have diminishing returns on health and thus having a stroke. For example decreasing an increasing rate. Or maybe (totally made up) exercising a little bit is really good for you, exercising a moderate amount is really bad for you, and exercising a lot is really good for you again. We might have a U shaped curve like (x-4)^2+4 for example when plotting health to exercise. I would image that having a curved, non-linear independent variable in a logistic regression like the one I described above would cause problems. Am I right? And if so, what kinds of things can you do to still include the variable in the logistic regression. Perhaps a transformation? Is that all?


title: R lme4: how to apply binomial GLM to percentages rather than yes-no counts

text: I have a repeated-measures experiment where the dependent variable is a percentage, and I have multiple factors as independent variables. I'd like to use `glmer` from the R package `lme4` to treat it as a logistic regression problem (by specifying `family=binomial`) since it seems to accommodate this setup directly. My data looks like this:                > head(data.xvsy)        foldnum      featureset noisered pooldur dpoolmode       auc      1       0         mfcc-ms      nr0       1      mean 0.6760438      2       1         mfcc-ms      nr0       1      mean 0.6739482      3       0    melspec-maxp    nr075       1       max 0.8141421      4       1    melspec-maxp    nr075       1       max 0.7822994      5       0 chrmpeak-tpor1d    nr075       1       max 0.6547476      6       1 chrmpeak-tpor1d    nr075       1       max 0.6699825      and here's the R command that I was hoping would be appropriate:                glmer(auc~1+featureset*noisered*pooldur*dpoolmode+(1|foldnum), data.xvsy, family=binomial)      The problem with this is that the command complains about my dependent variable not being integers:               In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!      and the analysis of this (pilot) data gives weird answers as a result. I understand why the `binomial` family expects integers (yes-no counts), but it seems it should be OK to regress percentage data directly. How to do this?


title: Sensitivity and specificity differ from my text box?

text: How is the sensitivity and specificity defined in the following confusion matrix from the caret package? According to my textbook ('Introduction to Data Mining' by Tan, Steinbach and Kumar) the sensitivity and specificity would be equal to the 'Pos Pred Value' and 'Neg Pred Value' found below: Moreover, which definition is used for making ROC-curves? I use the pROC package.               Confusion Matrix and Statistics                     y_BBR_test_est                 Mersalg Stilstand       Mersalg       432       309       Stilstand    2906       793                         Accuracy : 0.2759                                95% CI : (0.2628, 0.2893)         No Information Rate : 0.7518                   P-Value [Acc > NIR] : 1                                           Kappa : -0.0844               Mcnemar's Test P-Value : <2e-16                                Sensitivity : 0.1294                           Specificity : 0.7196                        Pos Pred Value : 0.5830                        Neg Pred Value : 0.2144                            Prevalence : 0.7518                        Detection Rate : 0.0973                  Detection Prevalence : 0.1669                     Balanced Accuracy : 0.4245                           'Positive' Class : Mersalg


title: What to do with missing values?

text: I have missing values for some of the variables in my data. I am using pooled OLS and have 144 observations. I have missing values for three of the variables. Less than 10% of the data for each variable is missing. What should I do with these missing values? Should I use multiple imputation or should I just omit the observations with missing values? Can I bias thw results if I use multiple imputation? * * * Thanks for your answer. I am using country data for the years 2010-2012, and the missing data is on trade openness and GDP per capita. Data on GDP per capita is missing for Syria and Cuba. Data on trade openness measured as (export+import)/GDP is also missing for Syria along with other countries, like Myanmar, Iran, Libya and for one year for Hungary, Cambodia, and Malawi. I can see that the missingness of data for Syria is not random, because of the war etc.The same holds for Libya and Iran. However, for trade openness there is also missing values for Hungary and Cambodia which I think is random (cant see other reasons why is should not be random).   These variables are not of interest for my regression, they are used as control variables.


title: Afex does not give me all p-values of all parameters. What should I do?

text: I have installed the package afex with this formula:               x <- mixed(duration ~ (1|item) + (1+color|speaker) + group*color*sex, data=data1.frame,REML=FALSE,type=2,method=c("KR"))      The formula gave me this:               Fitting 10 (g)lmer() models:     [..........]     Obtaining 7 p-values:     [Error in list2env(data) : first argument must be a named list     In addition: Warning message:     In mixed(duration ~ (1 | item) + (1 + color | speaker) + group *  :       Implementation of Type 2 method not unproblematic.       Check documentation or use car::Anova (Wald tests).      > x                              Effect stat ndf    ddf F.scaling p.value     1               group 1.61   2  35.74      0.98     .21     2               color 0.04   2  17.04      1.00     .96     3                 sex 9.71   1  66.83      1.00    .003     4         group:color 0.20   4  36.56      0.97     .93     5           group:sex 3.30   2 170.83      1.00     .04     6           color:sex 0.31   2  63.48      0.98     .73     7     group:color:sex 1.06   4 193.57      0.99     .38      Could I ask how I can get p values for all parameters?


title: What do the tests of model effects and parameter estimates really tell (when an interaction is defined)?

text: A couple of times in LMM or GEE (with SPSS, though I doubt that matters – and might occur in other analyses as well, but these are the ones with which I have seen it) I have seen something that seems contradictory: a continuous variable $X$ has completely different test values (Wald's or $F$s) and $p$s under tests of model effects and parameter estimates. For example, I have a GEE here with $p = .003$ for $X$ in tests of model effects, but $p = .728$ in parameter estimates: ![enter image description here](http://i.stack.imgur.com/MjU0L.png) I understand that these two $p$s for a categorical variable $Y$ can be different, as the tests of model effects check if any of the categories is different from zero, and if the reference category happens to be such that the estimate for category 1 is smaller and for category 2 bigger than the reference category, the estimates for those categories may not be different enough from the reference category (looking like I had no significant differences in parameter estimates) even if the difference between 1 and 2 would (which still shows in tests of model effects). This is clear. But I'm looking at a continuous variable now, with only one $b$ in parameter estimates that, if there were no interaction terms, has identical test values and $p$s with tests of model effects. However, when an interaction between $X$ and a categorical variable $Y$ is defined, the values may change wildly, and I don't know why. I assume it happens only in interactions because the parameter is estimated when the other parameters are fixed to zero (continuous) or to the last category (categorical), and when an interaction $XY$ is defined, the estimates for categories of $Y$ affects the estimation of $X$ (as can seen by the fact that the estimates for $X$ change when the reference category for $Y$ is changed). If I had only $Y$ in the model, I could change the model so that the most extreme category is the reference, and I would see that the tests of model effects would be reasonably similar to tests of parameter estimates. But interactions seem to work different somehow. When I change the reference category for $Y$, the tests for parameter estimates for $Y$, $X$, and $XY$ change (for $X, p$ changes from $.103$ to mentioned $.728)$, but in the tests of model effects, only the effect for $X$ changes, and by a much lesser amount than in parameter estimates (between $.001$ and $.014)$.   1. Why does the test for model effect for $X$ change when the reference category is changed? If the tests of model effects indicates whether the variable as a whole is different from zero, the reference category shouldn't affect it (like it does not affect tests for model effects for $Y$ and $XY$).   2. And why are the test values and significances so different? How does it make sense that there is an effect for $X$ in the model, but when I try to find what it is, the estimate for the parameter is always so small that it is not different from zero? At the same time, the model effects for $Y$ and $XY$ seem much more consistent when compared to parameter estimates. What am I missing?


title: Is adjusted R-squared appropriate to compare models with different response variables?

text: I heard that adjusted R-squared from two model are comparable only if two models use the same response variables. So if the response variable in one model is Y and the other one is log(Y), then I should not use adjusted R-squared; rather I should use AIC, BIC, etc instead. Is this true?


title: distortion function for k-means algorithm

text: I was reading Andrew Ng's ML lecture notes on K-mean clustering, in which the distortion function is defined as follow $$J(c,\mu) = \sum^m_{i=1} || x^{(i)} - \mu_{c^{(i)}}||^2$$ I am puzzled about the $L_2$ norm, since $|| x^{(i)} - \mu_{c^{(i)}}||^2 $ would imply $\sum^m_{i=1} (x^{(i)} - \mu_{c^{(i)}})^2$ and this means that there would be two summations $\sum_{i=1}^m$ in the entire expression. I am sensing that I have misunderstood something crucial here. Please point out the error. Thanks. UPDATE: the problem has a given a training set $\\{x^{(1)}, ..., x^{(m)}\\}$, where $x^{(i)} \in \mathbb{R}^n$ and the cluster centroids are $\mu_1, \mu_2,...\mu_k \in \mathbb{R}^n$


title: Smooth expectations outside the exponential family

text: At page 85-86 of Young and Smith "Essentials of Statistical Inference" there is an interesting result. If $X$ is a r.v. distributed according to the exponential family and $\phi(x)$ is a bounded (but not necessarily differentiable) function of x, then: $$ E(\phi(x)) = \int \phi(x) c(\theta) h(x) \text{exp} \big(\sum_{i = 1}^k \theta_i \tau_i(x)\big )dx $$ is differentiable of all orders w.r.t. $\theta$. My interpretation is that taking expectation under the exponential family smooths the function $\phi(x)$. I was wondering if similar "smoothing'' effects take place also outside the exponential family.


title: With a small sample from a normal distribution, do you simulate using a t distribution?

text: I want to simulate temperature data for some "what-if" calculations. The problem is, I only have a time series of 10 actual temperature data values. I want to use temperature as an input to the simulation, so I need a way to generate a large number of temperature values that are consistent with the original 10 values. It's probably ok to assume that they came from a normal distribution, but I don't know the mean or the variance. I have no way to prove it, but I doubt the 10 values do a good job of representing the full temperature range. If I use the `sample` function for the simulation, as shown below, I only get the original values back. That just doesn't look right. If I use the `rnorm` function, I know that I don't know the variance, so I don't think that is right either. So, I'm left with the `rt` function (t-distribution). Below is a mock up of the problem.               ori <- rnorm(n=10, mean=65, sd=5) #original 10 data points          num.sam <- sample(x=ori, size=100, replace = TRUE) #simulation using sample     num.tdis <- mean(ori) + (rt(n=100, df=10) * sd(ori)) #simulation using a t distribution          hist(ori, breaks=40:90)      hist(num.sam, breaks=40:90)      hist(num.tdis, breaks=40:90)       My questions are,   1. When I only have data (mean and variance unknown), and it is reasonable to assume that the data came from a normal distribution, is it ok to generate data for a simulation using a t-distribution?   2. For this type of situation, the only time I would use `rnorm` for the simulation is if I knew the variance (not a variance estimated from the data), right?   3. If a t-distribution simulation is ok for these conditions, are there any conditions where it is better to just sample the data (for example 100 original data points, 200, etc)? Edit:   1. Since I used the original data to estimate the mean and variance, should the degrees of freedom in the third line of the code (for `rt(...)`) be reduced from 10 to 9? Or 8?


title: Normality test for anova mixed model in SPSS: I did it in two ways but the results are contrasting

text: I'm conducting an ANOVA mixed model 2 X 3 (group X condition). Right now, I'm checking for the assumptions of ANOVA, such as normality. I should check indipendently for group and condition, right? I did it in two ways: 1) By, at first, splitting the sample (selecting "Data/Split file..." and entering the group variable) and then running the "1 sample Kolmogorov-Smirnov Test" and choosing "Normal distribution". 2) Without splitting the sample, but selecting "Analyze/Descriptive Statistics/Explore..." and entering the group variable as a factor and clicking "Plots..." and checking "Normality plots with test". The 1) returns all variables absolutely not significant (0.80 < p < 0.90). The 2) ran both the Kolmogorov-Smirnov and Shapiro-Wilk tests, with the former resulting not significant for all variables (but all p values lower than or equal to 0.20), instead the latter in some cases considerably significant (e.g. 0.012). I don't know which one should I trust, can you help me? Thanks!


title: How to handle strings in numeric data columns in a dataset?

text: In my data set there are columns like               Percent of Heart Attack Patients Given Aspirin At Arrival     ---------------------------------------------------------     Not Available     99%     Too few cases     100%     66%     Not Available      What is the generic process to handle such attributes which are numeric, but contain text like these? I'm just looking for what would be a standard acceptable approach that I can explain in my report - but if its a tool specific answer - I have rapidminer, R and IBM SPSS Modeler to work with.


title: Building Ensemble Models Practical Question

text: I understand the general idea of stacking several models and combining them to create a model that might perform better than any of the individual models that it is composed of. My question is regarding the actual steps needed to perform this. My thinking is that you would use the predictions of the individual models as features and use logistic regression (for a classification problem) to train them to get your ensemble model. The issue that I don't understand is what data to use. If you split your data into a training set and a test set and train your individual models on the training set and use the test set predictions as features for the ensemble model then that would mean that you must train your ensemble model on the test set. If this is the case then how would you test your ensemble model? What data would you use? Should you split your data into 3 sets - one to train the individual models, one to test the individual models and train the ensemble model, and finally one to test the ensemble model? I guess this ensemble concepts is confusing me as it relates to training sets and test sets. Thanks


title: Book recommendation for regression techniques

text: I'm taking my first course that contains regression techniques, analysis of variances and time series. What would you recommend if I would like to lean the mathematical theory behind those methods? For example, my lecture notes says that in least square method one can show that the estimation of regression coefficients has the form $$s.e.(\hat{\beta_0})=\sqrt{\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{n-2}\left (\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n (x-\bar{x})^2}\right )}$$ but I would like to learn proofs of such things.


title: Have several low order time series. How do I generate the best prediction?

text: I have a bunch of data in times series format, December 2013 through March 2014. I have this data for many different factors. Is there a way to take into account all 15000 or so observations in generating a prediction for each one for April 2014 even though each specific factor only has 4 data points?


title: Example of PCA where low variation PC is "useful"

text: Normally in PC the first few PC's are used and the low variation PC's are dropped -as they do not explain much of the variation in the data. However, are there examples where the low variation PC's are useful (ie. have use in the context of the data.. has an intuitive explanation and should not be thrown away)? Thanks!


title: Treatment of triple seasonal data

text: I have a data set of electricity spot prices, which contains three kinds of seasonality: one within 24 hours, one within a week and one within a year. I want to use an R package (`tsDyn`) which can't cope with seasonality, so first I would like to remove all three seasonalities, then adapt a model to the deseasonalized data, perform a forecast and then add the seasonalities, if it is possible, in order to transform my forecasts to reasonable form. Is this approach sensible and possible? And if yes, how could I accomplish this triple deseasonalization and then undo it within R? In the case of a simple one lag differencing I would just undo the seasonal differencing with 'cumsum()', but is something like this applicable for my data set?


title: How to compare forecasting methods: based on ARIMA and curve fitting?

text: I'm making a project connected with identifying the dynamics of sales. My database concerns 26 weeks (so equally in 26 time-series observations) after launching the product. I want to make forecast based on S-curve for clusters of time-series. The main aim was to compare two methods of forecasting: => based on parameters of logistic curve => based on ARIMA However, I do not know how to compare these two methods, because ARIMA bases on forecast differences and S-curve forecasts culminative sales. I tried with dm.test but S-curve has higher values of errors due to scale. So how can I compare which forecast method is better? I would be grateful for any explanation


title: How to train input feature weights for a path optimization problem?

text: I am working on a general path optimization problem. As many of you know, once the weights on all nodes are determined, one may solve this problem in many different means. Unfortunately, my raw input weights are often too noisy to be used directly. Therefore, I have to estimate the weights myself. I have been given some annotated data, where the optimized paths are manually drawn. I wonder how to effectively use these annotated data to improve my existing system. In particular, I have $N$ different denoising methods, each gives me an estimated node weights, or equivalently I have $N$ features for each node. Is it possible for me to train weights for each feature, such that my system works better in terms closer to the annotated data? By the way, I've already had some tools to evaluate the similarity between two paths. However, the direct use of these similarity scores as feature weights does not help, and the corresponding system is worse than my single best system. My guess says if I can find a reasonable set of weights for each feature, it is possible to improve my system. The difficulty of the problem lies here: if we have been given the clean weights on each node, then we could construct a direct function relationship between inputs and outputs, but now we only have the annotated optimized paths, which just reflects on these clean data to some extent. Can anyone point me some materials related to this problem? Many thanks.


title: is this interpretation of the p-value legit?

text: The interpretation of the p-value is a very difficult thing because what is legit and what is not is very close together. The p-value is the probability of the resulting test statistic if the null hypothesis $H_0$ is true; that means $\Pr(\mathrm{data}|H_0)$. A lot of people make the mistake of confusing this with $\Pr(H_0|\mathrm{data})$, which is wrong. Now to my question: I stumbled upon an online independent t-test and the result was a p-value of 0.20. The given interpretation was the following: "There is a 79.76% chance the proportions are different." My question: is this interpretation legit? Or is it the same confusion as mentioned above (e.g. it describes the probability of $H_0$). Thanks in advance!


title: Filter algorithm for a system

text: I have a system with the following structure: $X_{t+1} = X_t + E_{t+1}$ $E_{t+1} \sim N(0, \Sigma)$ $Y_{t+1} = f(X_{t+1})$ $Y_{t+1} \sim {\rm Uniform}(a_{t+1}, b_{t+1})$ So $X$ is a vector of hidden states. $E$ is an error. $Y_{t+1}$ is determined by $X_{t+1}$ via a deterministic non-linear function $f(.)$. $a_{t+1}$ and $b_{t+1}$ are observations at time $t+1$ and known. Given a list of observations $(a_1, b_1)$, ..., $(a_{t},b_t)$, I want to estimate $P(X_t|(a_1,b_1,a_2,b_2,...,a_t,b_t))$. Is there an efficient way to do this in real-time?


title: Polynomial distributed lag glms

text: I am about to use this input for my data but i cant find under which packages does it work. I searched almost everywhere. m<-pdlglm(yy~pdl(zz,4,2)+sin(zz)+pdl(xx,4,4))


title: Is dynamically increasing the sample size okay, if stated a priori?

text: I'm about to do a study about the merits of one stimulus compared to another with a within-subject design. I have a permutation scheme that is designed to reduce order effects of some parts of the study (task type order, stimulus order, task set order). The permutation scheme dictates that the sample size be divisible by 8. To determine the sample size I'd have to either take a wild guess (being a good tradition in my field) or calculate the sample size for my desired power. The problem is now that I don't have the slightest clue how large an effect size I'm going to observe (also a good tradition in my field). That means power calculation is a bit difficult. Taking a wild guess, on the other hand, might also be bad because I can either come out with too low a sample size or pay too much money to my participants and spent too much time in the lab. Is it okay to state upfront that I add participants in batches of 8 persons until I leave a corridor of two p-values? E.g. 0,05 < p < 0,30? Or how else, would you recommend, should I proceed?


title: How to identify outliers in server uptime performance data?

text: I have a python script that creates a list of lists of server uptime and performance data, where each sub-list (or 'row') contains a particular cluster's stats. For example, nicely formatted it looks something like this:               -------  -------------  ------------  ----------  -------------------     Cluster  %Availability  Requests/Sec  Errors/Sec  %Memory_Utilization     -------  -------------  ------------  ----------  -------------------     ams-a    98.099          1012         678          91     bos-a    98.099          1111         12           91     bos-b    55.123          1513         576          22     lax-a    99.110          988          10           89     pdx-a    98.123          1121         11           90     ord-b    75.005          1301         123          100     sjc-a    99.020          1000         10           88     ...(so on)...      So in list form, it might look like:               [[ams-a,98.099,1012,678,91],[bos-a,98.099,1111,12,91],...]      ### My question:   * What's the best way to determine the outliers in each column? Or are outliers not necessarily the best way to attack the problem of finding 'badness'?  In the data above, I'd definitely want to know about bos-b and ord-b, as well as ams-a since it's error rate is so high, but the others can be discarded. Depending on the column, since higher is not necessarily worse, nor is lower, I'm trying to figure out the most efficient way to do this. Seems like numpy gets mentioned a lot for this sort of stuff, but not sure where to even start with it (sadly, I'm more sysadmin than statistician...). When I asked over at Stack Overflow, someone mentioned using numpy's scoreatpercentile function and throw out anything over 99th percentile - does that seem like a good idea? (Cross-posted from stackoverflow, here: http://stackoverflow.com/questions/4606288)


title: Which multiple comparisons correction should I use?

text: I'm studying structure-function relationships in neurons and mainly focusing on correlations between dendritic length and excitability. My question is: **When searching for correlations, which multiple comparison correction should I use to classify which of the correlations I find are significant?** I have measured dendrite length (apical and basal) which provided me with 3 values: apical dendrite length, basal dendrite length and total dendrite length. I've also measured 5 different electrophysiological paramaters within the same cells. My research question is: Are the electrophysiological parameters related (correlated) with the dendritic length? This yields 3×5 comparisons (3 length parameters with 5 electrophysiological). I have found some significant correlations, as expected. Now I want to find out whether the significance of these correlations must be corrected, and if so, which method I must use. Bonferoni correction seems a little too conservative, because correlations between structure and function of a neuron seems logical, and it has already been established that they exist.


title: Improving accuracy of a binary classification when the target is unbalanced

text: I am working on the BRFSS dataset with the goal of predicting Diabetes. The dataset has 500,000 rows and 405 columns. It is a 0/1 classification problem, the ratio of 0 to 1 is 90:10. I tried using decision trees, logistic regression an ensemble of decision trees and logistic regression and my misclassification rate is almost 14% in all of these methods.   * What should I do to increase the accuracy?  I saw an earlier post which says subsampling or assigning different weights helps. But I am not sure about the ratio.   * What would be the best ratio to start off with?    * I am working using SAS. Is there a way to do subsampling in SAS?    * I am also interested in trying out the weighted approach. Is there a way to implement this in SAS? EDIT (28 Apr 2011) I tried subsampling and my misclassification rate goes up from 14% to 23%. The ratio I used was 50:50 for classes 0 and 1. The original ratio in the data was 90:10 and using the data as it is gave 14% error. So I believe subsampling doesn't work for my data. Would you suggest any other way to improve accuracy?


title: Can slopes in linear regressions be used as the independent or dependent variables in other regression models?

text: I have 100 patients and each patient have 10 longitudinal serum creatinine measurements. The estimated glomerular filtration rates (eGFR) were calculated from a MDRD formula comprising gender, age and serum creatinine. eGFR is the dependent variable and time is the independent variable in linear regression for each patient.   1. Do linear regressions violate the "independent X's" assumption and linear mixed models should be used instead?   2. Can eGFR slopes (which are estimates with uncertainties instead of measured numbers) estimated from each patient (in linear regressions for each patient or in random effects mixed models [how to estimate slopes for each individual patient in mixed models?]) be used as the independent or dependent variables in other regression models? Thank you.


title: Free public interest data hosting?

text: I have hourly and daily temperature reports for many stations at http://data.barrycarter.info/ I encourage people to download it, but, at 6.6G, it uses up a lot of bandwidth. Is there a service that hosts "public interest" data for free? I know about http://aws.amazon.com/publicdatasets, but you need an Amazon EC2 account to access that data.


title: If the empirical distribution of a sample is same as the true distribution, how shall the deviation be estimated?

text: Suppose we can always generate a sample $\\{X_1, \dots, X_n\\}$ of size $n$ for a discrete distribution, such that its empirical pmf is the same as the true pmf. One of the reason that we can have such assurance is when $X_i$ is a transform of a uniform discrete random variable $Y_i$, and we can generate a sample of $Y_i$ by deterministically enumerating the members of its range. for example if the common distribution is a Bernoulli distribution with parameter $p$. If we can always ensure that $\sum_i X_i /n = p$, shall the deviation be estimated as ${\frac{1}{n} \sum_i (X_i - \bar{X}}^2)^{-1/2}$ or ${\frac{1}{n-1} \sum_i (X_i - \bar{X}}^2)^{-1/2}$. Normally we would prefer the latter because it is unbiased while the former isn't. But here we know exactly that $\sum_i X_i /n = p$. Will the former be better than the latter? Thanks and regards!


title: Spatial interpolation models: deterministic vs statistical

text: I am applying diferent methods to interpolate continuous spatial surfaces (kriging, splines, glm,etc). Most of the studies that have enough detail for me to follow usually focus on one specific method. Is there any study that I might have missed where deterministic and statistical models are compared? I understand that the selection of the best model is made by the quantification of error and statistical techniques, what I am looking for are studies/examples that for a certain type of dataset(static, time-series, climate, topographic) show that one model is more adapted than other and why.


title: Different types of error rates when there are multiple comparisons

text: I have seen 5 different type of error rates. They include:   * **Componentwise Error Rate**   * **Experimentwise Error Rate**   * **False Discovery Rate**   * **Strong Familywise Error Rate**   * **Simultaneous Confidence Intervals** Suppose $H_0 = H_{01} \cap \cdots \cap H_{0K}$ Here are some questions that I have:   1. So basically if the componentwise error rate is $\alpha$ then the probability of rejecting a single hypothesis in a single test is $\alpha$?    2. If the experimentwise error rate is $\alpha$ then the probability of rejecting any of the $H_{0i}$ when all of the $H_{0i}$ are true is $\alpha$?   3. What is the difference between the FDR and the experimentwise error rate? Doesn't the FDR control the experimentwise error rate?   4. Is the strong familywise error rate basically the strictest error rate? So the FDR allows for some slack if we have more correct rejections but the strong familywise error rate does not?   5. Simultaneous confidence intervals must cover their true parameter with probability $1-\alpha$. But a single confidence interval will cover its true parameter with probability greater than $1-\alpha$?


title: Uncertainty from linear fit on additional data

text: Let's say I have 5 known data points with coordinates $X$ : Area under curve $Y$ : Activity The 5 points have individual error ($\Delta X_{i}$,$\Delta Y_{i}$) on both $X$ and $Y$ and I know that the best fit of these 5 points is a linear fit. After plotting my 5 original points, I have another one that I only know the Area under curve and not the activity. With the fit equation of my 5 original points, I can associate an Activity to my 6th point. Is there a way to provide a uncertainty to my 6th point ? I have some ideas about Max/Min slopes but I don't know if it's the best way in this situation... And I have to use only Excel at the request of someone. Thank you !!


title: Optimal bin width for two dimensional histogram

text: There are lots of rules for selecting an optimal bin width in a 1D histogram (see for example) I'm looking for a rule that applies the selection of optimal equal-bin widths on _two-dimensional_ histograms. Is there such a rule? Perhaps one of the well-known rules for 1D histograms can be easily adapted, if so, could you give some minimal details on how to do so?


title: Multiple regression with class probabilities

text: I'm given a list of class probabilities $(p_1, p_2, \ldots, p_c)$ (which sum to 1) for each sample. Now, I'd like to replicate this with a model (supervised learning). How can I approach this? I could use a separate regression model for each of $p_1, \dots, p_{c-1}$. But the probabilities might exceed 1 and also the models are independent, which both loses information. I could use a classification method where I try to predict probabilities for classes. But then the input for each sample is a single class (the largest class probability). Here I'm losing the details of the given probabilities. The goal is actually to closely match the given class probabilities. Anything more direct and appropriate?


title: Half of the mean as biased estimate - Stein's paradox

text: i read this paper with interest yesterday http://statweb.stanford.edu/~ckirby/brad/other/Article1977.pdf and stumbled upon the statement that, if and only if, the population mean is close to zero, than the risk (MSE) of using half of the sample mean as an estimate for the population mean is lower than that of using the mean. I tried to model this out (unfortunately in Crystal Ball / Excel, did not use R for this) but could not replicate this result. in my examples, the risk (as per MSE) of the half- of-the-mean estimate became closer to the mean estimate, but never became lower. I might, clearly, have misunderstood the concept, but would be very interested if someone can show / further explain this to me, Many thanks in advance!


title: Example where a simple correlation coefficient has a sign opposite to that of the corresponding partial correlation coefficient

text: > Give some examples where a simple correlation coefficient has a sign > opposite to that of the corresponding partial correlation coefficient and > comment on it. It is a question from an examination paper.Please help.


title: Standardized solution for censored regression

text: The censReg package in R provides identical results to Mplus for censored regression. Mplus also provides a fully standardized solution. How can I scale the censReg coefficients so they are standardized, (beta weights)? Does it make sense to standardize marginal effects in censored outcome models, (see the margEff.censReg function)?


title: Differential Equations as Generative Models

text: I wonder, if we can say stochastic differential equations _are_ generative models. I usually think about the Kalman filtering for example, we fix a discrete-time evolution equation of a certain object, which is a differential equation in continuous domain, and then use this model as a generative model to derive inference algorithms in a Bayesian setting. Then, can we say that all differential equations involving noise terms are generative models of some underlying phenomenon? PS (December 7, 2012): Although it is discussed in the comments, I push my luck once more by adding a concrete example. Consider a stochastic difference equation, $$x_{t+1} = x_t + w_t$$ where $w_t \sim \mathcal{N}(0,1)$. Then, we can write in the probabilistic terms, $$p(x_{t+1}|x_t) = \mathcal{N}(0,x_t)$$ and this is a generative model. Of course, we can write another equation, such as an observation model, and this picture completely define a discrete-time stochastic dynamical system as well as a generative model. Then I reask my question: is there a specific reference that promotes this relationship? Or is it quite trivial and obvious for researchers from both sides? Thanks.


title: SVM retrain on whole dataset for final model --> overfitting?

text: i am training a SVM (RBF kernel) with a dataset of ~1500 samples (balanced) using fminsearch on the CV error for parameter optimization (C and s). After i found the "best" parameters (local optima possible) i am retraining my model on the whole dataset to derive a "final" model. Is this a wise thing to do? Would the final model be proned to overfitting? I experience worse performance on unseen data which might be OK as the CV during my optimization approach produces a somewhat optimistic estimate on the test error. I think this adresses a pretty general problem but i could not find proper reasoning yet... Would it be a reasonable alternative to use just one of the models from the best performing crossvalidation? I assumed that once the parameters are fixed the model will not suffer from overfitting no matter how many more samples i use for training? Thanks, Pir


title: Small sample MLE

text: Is there a go-to reference for the effect of small samples on MLE estimates? Specifically the effect on the parameter distribution. i.e. should we switch to something like a _t_ -distribution with fewer degrees of freedom rather then a normal for the distribution of the parameters?


title: prior distribution to the binomial distribution probability distributions urn model

text: I have an infinite population with unknown mean of successes and failures. I'm drawing 400 times from the population and get 400 successes. Now I want to generate random estimates for the true mean of the population from which I have drawn the 400 successes. Can anyone tell me which probability distribution I have to use, respectively how the function to use in R would look like? I thought that               rbeta(1,400+1,400-400+1)      might be the right function, but even if I perform this one 10 000 000 times I never get the result 1 (and the same accounts for `rbeta(10000000,0+1,400-0+1)` and the result 0). So I ask myself why isn't it possible to have a population which consists only of successes or only of failures? Many thanks in advance!