title: Calculation of canonical link function in GLM

text: I thought that the canonical link function $g(\cdot)$ comes from the natural parameter of exponential family. Say, consider the family $$ f(y,\theta,\psi)=\exp\left\\{\frac{y\theta-b(\theta)}{a(\psi)}-c(y,\psi)\right\\} $$ then $\theta=\theta(\mu)$ is the canonical link function. Take Bernoulli distribution as an example, we have $$ P(Y=y)=\mu^{y}(1-\mu)^{1-y}=\exp\left\\{y\log\frac{\mu}{1-\mu}+\log{(1-\mu)}\right\\} $$ So, the canonical link function $$g(\mu)=\log\frac{\mu}{1-\mu}$$ But when I see this slide, it claims that $$ g'(\mu)=\frac{1}{V(\mu)} $$ Though it can be easily verified for this particular distribution (and some other distributions, like Poisson distribution), I can't see the equivalence for the general case. Can anyone give hints? Thank you~


title: Individual level prediction of a person’s probability of voting without their vote history

text: Is it possible to create individual level predictions of a voter's probability of voting when you do not know their vote history? In the data set provided in my homework assignment, I am given data on a number of people. Some have a known voting history and some have missing values within the variable (these are the people the people without known voting histories). The voting history is labeled `voted` with binary values (0,1) and the unknown vote histories show up as missing values (in stata: `.`). I am also provided with other variables in the data set such as sex, age, income level, labor force, and race. I have tried replacing the missing values and running a logistic regression against all of the other independent variables such as sex, race, age, income, level and it has not worked. _UPDATE: What I have tried thus far:_   I have also tried to create a new variable with just the missing variables in an attempt to isolate the people who have no vote history to make the problem easier but that also has not worked.


title: How to prove that test 1 produces more results than test 2?

text: If I have 2 tests applied on a set of data, both tests either produce a result ( assume = 1) or no result ( NA , assume = 0 ). Is there a way to prove that test 1 produce more results than test 2 with significant results? For example the results below shows that test 1 is better than test 2 but how to prove the significance part ? Is there other test than using t-test?               Test1 , test 2     1, 0      1, 0      0, 0,      1, 1     0 , 0      0 ,1     1, ,1     1, 0      1, 0       Thanks


title: Coin flipping, decision processes and value of information

text: Imagine the following setup: You have 2 coins, coin A which is _guaranteed_ to be fair, and coin B which may or may not be fair. You are asked to do 100 coin flips, and your objective is to _maximize the number of heads_. Your prior information about coin B is that it was flipped 3 times and yielded 1 head. If your decision rule was simply based on comparing the expected probability of heads of the 2 coins, you would flip coin A 100 times and be done with it. This is true even when using reasonable Bayesian estimations (posterior means) of the probabilities, since you have no reason to believe that coin B yields more heads. However, what if coin B is actually biased in favor of heads? Surely the "potential heads" you give up by flipping coin B a couple of times (and therefore gaining information about its statistical properties) would be valuable in some sense and therefore would factor into your decision. How can this "value of information" be described mathematically? **Question:** How do you construct an optimal decision rule mathematically in this scenario?


title: Mann-Whitney U Test Analysis Help

text: Good day! I'm totally new to SPSS (and the test as well, but after some research I've finally demystified a bit about Mann-Whitney and how to do the test in SPSS. The only problem I'm currently facing is the analysis) And so in this snippet, I want to discover how the researcher came up with this analysis: ![Valid XHTML](https://fbcdn-sphotos-g-a.akamaihd.net/hphotos-ak- ash3/t1.0-9/q71/s720x720/10268622_10202139109587436_5776841905445543571_n.jpg). How was he able to come up with a conclusion that there is no significant difference between the examination marks? What specific data in the table made him conclude that? What variables did he look at? If I were to use a Mann- Whitney test and have to interpret the results, where will I look at? What are the factors do I have to consider? Any help would be very much appreciated. :) Also, what's the purpose of the Z and Wilcoxon W? Thank you very much!


title: Binary and ordinal scale

text: I know my response variable is ordinal(5 levels). I decide to ignore all the information in the response variable and group it into two groups. I then fit a logistic regression model to it. I use this logistic regression model for classification for the two groups. Alternatively, I fit an ordinal logistic regression model, with this ordinal logistic regression model I use it for classification for my outcome but only for classifying objects in the predefined groups that I choose before fitting the binary logistic regression model. Which model will perform better when assessed with the $AUC$?.


title: Why change items which don't change Cronbach Alpha?

text: I am doing peer review on a paper submission for a journal, and I am somewhat out of my depth. The authors are developing a psychometric scale. They did a pilot study with 45 items. They measured Cronbach alpha and found that it is 0.89 overall. Then they found that six items have `Cronbach's alpha if item deleted` of 0.89 or very similar (between 0.89 and 0.893). They show a table with all six. They then write > A number of modifications were made according to what Table 2 suggested And go on to explain how the wording of each item was changed to make it clearer. Obviously, they didn't want to have items whose `Cronbach alpha if deleted` is the same as the overall Cronbach alpha. Why are they doing this? Is this a common technique in constructing psychometric scales? Are there guidelines for doing it? They also don't cite any literature about the method they followed to develop their scale. It does look sound overall (pilot study, exploratory factor analysis, then refining the questionnaire based on the results), but I'm not sure whether the decisions they make about details (such as the above change of six items) are justified. Is the knowledge how to make a scale so common that there is no need to cite sources?


title: Individuals standard deviations and/or standard errors for groups after implementing ANOVA?

text: I've seen couple of articles reporting individuals standard deviations and/or standard errors for groups even after implementing ANOVA. My understanding is that groups SE's should be based on experimental error mean square. Any comment?


title: Obtaining $p$-values for UMP $t$ tests for equivalence

text: I am trying to understand the uniformly most powerful (UMP) equivalence tests presented by Welleck (2010) with respect to $p$-values for the tests. In the case of the paired $t$ test example concerning measures before and after an intervention (pp 94–96) Wellek presents a UMP one-sample $t$ test of mean equivalence with _symmetric_ equivalence boundaries where $ -\theta_{1} = \theta_{2} = \varepsilon = 0.5$, for the negativist null hypothesis H$^{-}_{0}\text{: }\left|t\right| \ge \tilde{C}_{\alpha,\nu}(\varepsilon)$ and gives:   1. $n=23$ (sample size)   2. $t=\frac{0.16}{\frac{3.99}{\sqrt{23}}} = 0.1923$ (test statistic)   3. $\tilde{C}_{.05, 22}(0.5) = \sqrt{F^{-1}_{\nu_{\text{n}},\nu_{\text{d}},\psi^{2}}(0.05)} = 0.7595$ (critical value), where $F^{-1}$ is the quantile function (inverse CDF) of the noncentral _F_ distribution, $\nu_{\text{n}}= 1$, $\nu_{\text{d}}= n-1 = 22$, and $\psi^{2} = n\varepsilon^{2} = 5.75$. We would therefore reject H$^{-}_{0}$ in favor of H$^{-}_{\text{A}}$, and conclude that, for $\alpha=0.05$ and $\varepsilon = 0.5$, mean measurements before and after intervention were equivalent (we would also draw this conclusion in a combined test of H$^{-}_{0}$ and H$^{+}_{0}$, since we would reject the first, but not the second). What is the $p$-value of this test? My guess is that it is **[edited to remove the "$\mathbf{1-}$" confusion I had... forgot which tail I was in]** :   4. $p = \text{P}(\left|T\right| < \left|t\right|) = F_{\nu_{\text{n}},\nu_{\text{d}},\psi^{2}}(t^2) = 0.00882$, where $F$ is the CDF of the noncentral _F_ distribution. Is this the correct $p$-value for this specific example, and is this the way to approach $p$-values generally for _t_ tests for equivalence? In the non-symmetric equivalence interval case, Wellek describes an iterative process to numerically solve for the values $C_{1}$ and $C_{2}$, corresponding to the rejection regions H$^{-}_{0}\text{: } t \le C_{1}$ or $t \ge C_{2}$. Any suggestions for how to obtain the $p$-values for $t$ in this case? **Update:** would: $\psi = n\left|\theta_{1}\theta_{2}\right|$ in (4) above in the asymmetric case? **Update 2:** Bounty for explicitly describing the steps required to iteratively estimate $p$-values in the asymmetric case as hinted at in Horst's response, _including_ how to solve $\min_{i \in l,u}\left|T−t_{i}\right|=0$.    **References** Wellek, S. (2010). _Testing Statistical Hypotheses of Equivalence and Noninferiority_. Chapman and Hall/CRC Press, second edition.


title: Estimates and C.I. of percentiles for a survival function

text: In S-plus estimates of percentiles for a survival function can be obtained using the `qkaplanMeier` function (on the results of a call to kaplanMeier) like that:               kfit <-kaplanMeier(censor(TIME,STATUS)~1)     qkaplanMeier(kfit, c(.25, .5, .75))      How can I do this in R?. Those functions do not exist anymore. What if I also want the (asymptotic) confidence intervals for the percentiles? How can I get the mean for the survival time?


title: Chi-squared vs Fisher's Exact Test w/ 5x6 contingency table & some cells <5?

text: I don't have the minimum required data of 5 in each cell for a chi-squared test, but my contingency table is greater than 2x2, so I can't use a Fisher's exact. Would it help to use a chi-squared test set to either exact or Monte Carlo to prevent it from making the assumption of normal distribution? I'm looking at trauma patterns created by five different tools; six categories of trauma have been defined. I would like to see whether the patterns of trauma (by occurrence of each type) differ significantly from one another. N=67, but not all tools have all types of trauma.


title: Correlate 3 vectors that are sampled from 3 different distributions

text: I have three vectors with sampled values in MATLAB, each from a specific probability distribution. (The vectors contain time series data of the growth rate of retail prices of specific goods)               %pearsrnd(mu,sigma,skew,kurt,m,n)     commodity1=pearsrnd(0.005,0.085,0.237,7.899,1,600);     commodity2=pearsrnd(0.003,0.040,0.280,5.630,1,600);     commodity3=pearsrnd(0.006,0.139,-0.207,4.209,1,600);      However, I want to add two correlation to the data. Specifically, the correlation between commodity 1 and commodity 3 should be 0.135, the correlation between commodity 2 and 3 should be 0.115. How can I solve this most easily?


title: How to fetch all historically available hourly weather data?

text: Where can I obtain all hourly weather data available? Notes:   * Ideally, this would include the history of all data currently published at:      * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sflnd/DS.synop/     * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sflnd/DS.metar/     * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sfmar/DS.dbuoy/     * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sfmar/DS.ships/   * I've obtained 10 years of METAR data from wunderground.com like this:               curl -H 'Cookie: Prefs=|SHOWMETAR:1|;' -o data.txt 'http://www.wunderground.com/history/airport/KABQ/2005/05/02/DailyHistory.html?req_city=NA&req_state=NA&req_statename=NA&theprefset=SHOWMETAR&theprefvalue=1&format=1'       but am sure there's a better way.   * I'd love per-minute/etc data too, though I don't think anyone kept measurements better than hourly back then.


title: How to calculate a percentile of y for a given x given a series of (x, y)?

text: I have a scatter plot of (for example) height against age. How does one calculate for an individual point the percentile of the height for a given age? Suggestions in R would be most appreciated. Thanks!


title: How to compute term frequency and find clusters in a dataset composed of strings?

text: I am currently looking for some Information Retrieval techniques. I have a SQL database table containing strings. It has 1000 records, each being a random sentence I picked from random web sites. I need to get the _term frequency_ and represent each string into a vector. I also need to cluster the records, e.g. using k-means. Does anyone know what is the best way to do this? Are there any tools I can use? I am new to this and looking for a jump off point.


title: Is a subset $B= {(-\infty,t),t\in \mathbb{R} }$ an event of sample space $\mathbb{R}$?

text: The sample space $\omega$ is $\mathbb{R}$. Here is a subset $B= {(-\infty,t),t\in \mathbb{R} }$. So $B \in \beta(\mathbb{R})$. Is this subset an event? How to show that ? If it is an event, is it because it is Borel?


title: Crosstab factor scores generated by factor analysis

text: i just wanted to know if i could use factor scores and crosstab it with demographics (i.e. gender, age, etc.)? I have 69 likert-scale variables and run factor analysis on SPSS. It gave me 10 new variables (types of personality. i just wanted to see the demographics of each new variables. Thanks!!!


title: Issues with feature selection in matlab

text: I am trying to use sequentialfs to do some feature selection in matlab. I have huge dimensional data of 22215 features. When I tried to use sequentialfs with svm as classifier so that it selects the best subset of features, it just keeps on running, probably its because of the huge number of dimensions. However, weka does the same thing very quickly. Weka has wrapper filters and it does it so quickly. What sort of heuristic does it use? Even though I tried with 5000 features in weka since it was not taking 22215 features, it gave me results quickly with wrapper filters. What should I do with sequentialfs in matlab This is the command I am using in matlab               c = cvpartition(yS1,'k',12);     opts = statset('display','iter');     [fs, history] = sequentialfs(@SVM_class_fun, X, yS1,'cv', c, 'options', opts);          SVM_class_fun          function err = SVM_class_fun(xTrain, yTrain, xTest, yTest)       model = svmtrain(xTrain, yTrain);       err = sum(svmclassify(model, xTest) ~= yTest);     end      Here the dimension of X is 100x22215 where I have 100 examples each of dimension 22215.


title: How are zero values handled in lm()?

text: I am working on a linear regression with R and there are many 0 values in my predictor variables. How are these handled in R's `lm()` function? Should I remove this data for more accurate analysis? Any advice is appreciated. Thanks.


title: Marginal effect calculation in Stata

text: My estimated dynamic GMM model results are as follows. $$Y = -0.009 + 0.009PF + 0.015PS - 0.001(PF\times PS)$$ I want to calculate the marginal effect of $PF$ at $PS$ (min, mean, max), which I can express as follows: marginal effect of $PF = 0.009 - 0.001 PS$. I can calculate it manually, but it is time-consuming. Is there any command that calculates the marginal effect after dynamic GMM along with its significance in Stata 12 other than `margins`, or any way that `margins` can work? The `margins` command in Stata does not work after dynamic GMM.


title: Mahalanobis distance between two bivariate distributions with different covariances

text: The question is pretty much contained in the title. What is the Mahalanobis distance for two distributions of different covariance matrices? What I have found till now assumes the same covariance for both distributions, i.e., something of this sort: $$\Delta^T \Sigma^{-1} \Delta$$ What if I have two different $\Sigma$s? Note:- The problem is this: there are two bivariate distributions that have the same dimensions but that are rotated and translated with respect to each other (sorry I come from a pure mathematical background, not a statistics one). I need to measure their degree of overlap/distance. * _Update: *_ What might or might not be implicit in what I'm asking is that I need a distance between the means of the two distributions. I _know_ where the means are, but since the two distributions are rotated with respect to one another, I need to assign different weights to different orientations and therefore a simple Euclidean distance between the means does not work. Now, as I have understood it, the Mahalanobis distance cannot be used to measure this information if the distributions are differently shaped (apparently it works with two multivariate normal distributions of identical covariances, but not in the general case). Is there a good measure that encodes this wish to encode orientations with different weights?


title: evaluating effect of treatment with repeating subjects

text: We have a situation where multiple subjects who are either treated or controls are involved in multiple (different) activities A1..A20. Each subject might be involved in any subset of the activities and the same measurement is taken during the activity. We would like to match treated with controls based on age, gender, and activity. What is the best way to avoid the biased caused by looking at the same subject multiple times (each time under a different activity). Some of the options we are considering: a. aggregate all activities per subject and then match and compare. Problem: each subject is involved in a different subset of activities. b. same as a, but try to match subjects based on age, gender, and subset of activities for which data is collected for the subjects. c. match treated and control activities by age,gender,the activity and do not use a subject more than once (throwing away much of the data). d. same as c, and use subjects multiple times.


title: Performing CART using unlabeled data

text: I am involving in a project of customer churn prediction. I have 10000 customers data. All the data are unlabeled data. Now my question is can I perform CART analysis with unlabeled data?


title: Interpreting vecm result with 4 cointegrating equations

text: I have performed Johansen cointegration test and the result shows that I have more than one cointegration for all my models. Below is the cointegration test result for one of my models (I have 6). Prior to this, I have performed unit root test and found that all variables are integrated of order (1). I also perform lag length criteria test to determine the optimal lag length for my VAR model. From my understanding, I need to use VECM to determine the long run equilibrium relations (i.e which variables are statistically significant in the long run equilibrium). Since I have 4 cointegration, the results show the cointegration equations for the top four variables listed in the results. Therefore, how can I determine the relationship among the 4 variables/ how the top 4 variables related to each other? For example, the dependent variable in this model in lnCA_CV (balance in current account) and I would like to know how the interest rates and macroeconomic variable affect lnCA_CV. Below are the VECM result for this model. ![VECM](http://i.stack.imgur.com/1FD78.png) ![Error correction](http://i.stack.imgur.com/zPoCA.png) Thanks a lot!


title: Beginner question: finding correlations with percentage data

text: I am trying to work out whether mobile, desktop or tablet users are more likely to perform an action on a site. The data I have is structured as follows:               week     clicked      mobile    tablet    desktop     1        1            104       97        205     1        0            204       214       348     2        1            128       108       257     2        0            207       222       360      ... though I have about 25 weeks' worth and the numbers are larger! I can munge the data and restructure it if needed (e.g. into percentages) but that is the structure of the raw data. I am relatively new to statistics and trying to educate myself and ensure the approach I take is valid and meaningful. My null hypothesis is "the type of device used makes no difference to the outcome", but I'd like to investigate whether that's not the case. What would be a valid statistical test to use in this scenario? I'm working in R, can also use pandas. Any advice gratefully appreciated. UPDATE: The research I've done suggests that I should restructure this into rows for mobile, desktop and tablet, with the percentage conversion for each row, and then conduct a chi-squared test. Does that sound reasonable?


title: Variable Selection Methods in R

text: `regsubsets` and `stepAIC` are the two most common options for variable selection in R; they can be found in the `leaps` and `MASS` packages respectively. My questions are:   * Statistically; What are the pro's and con's of each method   * When should you choose one method over the other    * How can you compare the results of these two methods


title: How to improve standard visualizations of a univariate distribution?

text: I have a data set with about 15,000 labeled observations of a single continuous value. What is the best way to plot this type of data? I'm playing around with various histograms and density plots, but I can't seem to figure out the best way to plot this data set. Any suggestions? My data looks like this:                label     value       -------   -------         foo       1.2         bar       6.2         baz       0.2         qux       4.7         ...       ...        This data set contains 15,000 values, each with a unique label. I am looking as to how best to create a visualization of the distribution of the data and see outliers. Here are two candidate plots I've generated. Both simplified the data more than I would like. Are there any additional ways I could plot the data and somehow integrate the labels into this plot? ![density plot](http://i.stack.imgur.com/MNqcS.png) ![boxplot](http://i.stack.imgur.com/MXQ2y.png)


title: Find the equation from generalized linear model output

text: Another basic question. Let's say that I generate the probability of an outcome based on a certain factor and plot the curve of that outcome. Is there any way to extract the equation for that curve from R.               mod = glm(winner ~ our_bid, data=mydat, family=binomial(link="logit"))     summary(mod)          > summary(mod)          Call:     glm(formula = winner ~ our_bid, family = binomial(link = "logit"),          data = mydat)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -0.7443  -0.6083  -0.5329  -0.4702   2.3518            Coefficients:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -9.781e-01  2.836e-02  -34.49   <2e-16 ***     our_bid     -2.050e-03  7.576e-05  -27.07   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 42850  on 49971  degrees of freedom     Residual deviance: 42094  on 49970  degrees of freedom     AIC: 42098          Number of Fisher Scoring iterations: 4          all.x <- expand.grid(winner=unique(winner), our_bid=unique(our_bid))     all.x     won = subset(all.x, winner == 1)          y.hat.new <- predict(mod, newdata=won, type="response")     options(max.print=5000000)     y.hat.new          plot(our_bid<-000:1000,predict(mod,newdata=data.frame(our_bid<-c(000:1000)),type="response"))      ![enter image description here](http://i.stack.imgur.com/WsP47.png) So can I go from this probability curve to an equation in R. Something like:               Probability = -0.08*bid3 + 0.0364*bid2 - 0.0281*bid + 4E-14      Thanks for help!


title: Anomaly detection of web browsing sequences

text: Please consider that I'm quite new to machine learning. I need to create models based on browsing patterns of web users and find deviations from that model. I'm using web server access log files. For ex: suppose I'm having data for following two users for a period of month: user1 = [p1,p2,p3] user2 = [p1,p2] Then I upload one days log(preprocessed) and need to check whether there are any deviations. Please not user identification is done. Need advice on this?


title: Conditional distribution of a discrete random variable given a continuous one

text: I can´t solve a conditional distribution problem using the simple Bayes´ Theorem for PDFs. The thing is, I know that `X|Y ~ Poisson(Y)` and that `Y ~ Gamma(α,β)`. I´m being asked to find the PDF for Y|X but I´m complete clueless about how to do so. Can anyone point me in the right direction?


title: Alternative to Kolmogorov-Smirnov test when parameters are estimated from the data

text: I need to compare whether two distributions are similar when the values are scaled by the mean of each of the distribution. One limitation of ks-test as per http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm is "if location, scale, and shape parameters are estimated from the data, the critical region of the K-S test is no longer valid." Consider for example: `Data1` consist of 10000 numbers from uniform random distribution [0,1] with mean 0.5 `Data2` consist of 10000 numbers from uniform random distribution [0,10] with mean 5.001 If I compare Data1 with Data2/10 then ks-test gives that both the distribution are same; while comparing Data1/0.5 with Data2/5.001 gives that the distribution are different. Is there a way to check the similarity between the distributions in such cases? **Edit:** As the answer suggests I can use ks-test where the p-value is determined via permutation. My additional difficulty is that the data-points are integers: `Data1` consist of 10000 integers from uniform random distribution [0,10] with mean 5 `Data2` consist of 10000 integers from uniform random distribution [0,100] with mean 50.001 Is there a test to compare whether `Data1` and `Data2` are similar apart from the scale? Further, I do not know the actual scale and I am determining it from the data. These examples are just a proxy for my actual data, which are two experiment where 10000 people rated a movie on a scale [0,10], while in other case 10000 different people rated the same movie on a scale [0,100]. I want to check apart from the scale can one say that whether the distributions are same or not.


title: Poisson Process arrivals

text: This is a homework problem. Between 10 AM and 6 PM visitors arrive at the Tate Modern Gallery in accordance with a Poisson process at the rate of 6 per minute. Determine the probability that 10 visitors arrive between 1:00 PM and 1:02 PM, given that 800 visitors arrive between noon and 2:00 PM. My attempt: I know that the 10 arrivals (in permuted order) are identically and independently, $S_1,\ldots, S_{10} \sim \text{Uniform(0,120)}$ so $P(\text{Arrival between 1:00 PM and 1:02 PM}) = \displaystyle\int_{60}^{62} \dfrac{1}{120} dt = \dfrac{2}{120}$ So my guess is that $$P(\text{10 arrive between 1:00 PM and 1:02 PM}|N(120)=800) = \displaystyle{800 \choose 10}\left(\dfrac{2}{120}\right)^{10}\left(1-\dfrac{2}{120}\right)^{790}$$ Haha, I think I have figured it out as I typed out the question. Do you guys think this looks right?               > dbinom(10, size=800, prob=1/60)     [1] 0.0791953


title: Can we visualise regression results through diagram?

text: I read a journal article, in which the authors made a diagram based on regression results (beta weights). Pretty much looks like a path analysis: ![](http://i.stack.imgur.com/n1XAi.jpg)   I'm very sure that path analysis/SEM was not employed as can be read from their paper: > Both logistic regression and multiple regression analyses were used. > Logistic regression was used for the prediction of condom use because the > dependent variable was dichotomous. Prediction of the two intention > variables was accomplished using multiple regression. This has the advantage > over logistic regression of providing estimates of the effect size. Is that OK? I'm not trying to judge the quality of the paper. On the contrary, if that practice is acceptable, I would like to use it in my paper. I wanted to do path analysis, but I'm afraid my sample size (66 subjects) is too small for that analysis, so I used regression instead.


title: Estimate values below the limit of detection from a left-censored dataset

text: I have 32/180 determinations of a laboratory marker below the limit of detection (LOD) and would like to estimate them by maximum likelihood estimation. By using Minitab's reliability / survival analysis I was able to find a well-fitting model (lognormal) and to estimate mean, SD, and percentiles. How can I estimate the single missing values that are below the LOD? Is it OK to just use the percentiles for the first 32 values (up to the 17.77th percentile)?


title: Confidence intervals for Cohen's h effect size

text: When examining the difference between two proportions I typically use Cohen's h (i.e, the difference between two arcsin-transformed proportions) for the effect size. Does anyone know how I could calculate 95% confidence intervals for Cohen's h?


title: Is there an R package to check the proportional odds assumpion?

text: is there an R package/function to implement the parallel slope test given a clm (ordinal package) or polr (MASS) package?


title: EFA: Can I remove/drop variables with non significant loadings and re-run the EFA?

text: I am employing EFA to 56 items. However, there were cross-loadings occurred and, therefore decision to drop the items is made. The question: The rotated components matrix showed there were a few items with no-significant loadings in any of the components, so, should I remove/drop the non-loading items and re-run the EFA until all the items were loaded into the respective component? Thanks


title: Test for significance in a time series using R

text: I've got a time series data set that looks at population numbers on a yearly basis, and I would like to know if there is a significant difference between certain years, and if so- which years are statistically similar to each other. I'm using R, so I'm hoping to get some advice on R code to do this. Here is simplified example time series looking at a population by year:               Year<-c(2001,2002,2003,2004,2005,2006)     Pop<-c(1,4,7,9,20,21)     DF<-data.frame(Year,Pop)      So given the above data, ideally I would hope for some method that would calculate if there is a significant trend over time (p value) and which years are statistically different from each other (I'm guessing 2001-04 would be in one grouping and 2005-06 would be in a second).


title: Comparing two nominal-scale response distributions

text: I have data from an experiment that consist of nominal-scale responses that subjects made in two treatments (the treatments are between subjects, not repeated measures). I've provided my data in data frame and table form below. My goal is to compare the distributions of the responses in the two treatments - `t1` and `t2` \- to see if they are significantly different from one another. I am not interested in differences among levels of the response variable (e.g., between response `a` and `c`). I'm unsure how to approach this. My first thought was a Fisher exact test of the 8x2 contingency table (`tab` below), which yields: `p-value = 0.2102`. However, it seems that this tests whether the response variable is associated with the treatment variable as a whole, rather than testing for a difference between the treatments. Any suggestions for a better approach are most welcome.               # data frame with treatments and responses     df <- structure(list(id = 1:30, treatment = structure(c(2L, 1L, 2L,          1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,          1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("t1",          "t2"), class = "factor"), response = structure(c(2L, 1L, 1L,          1L, 1L, 1L, 6L, 8L, 1L, 1L, 4L, 1L, 1L, 1L, 5L, 1L, 7L, 1L, 8L,          1L, 4L, 1L, 4L, 2L, 1L, 1L, 1L, 3L, 1L, 1L), .Label = c("a",          "b", "c", "d", "e", "f", "g", "h"), class = "factor")), .Names = c("id",          "treatment", "response"), row.names = c(1L, 2L, 3L, 4L, 5L, 6L,          7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L,          20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L), class = "data.frame")          # 8x2 contingency table     tab <- structure(c(12L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 7L, 1L, 0L, 3L,          1L, 1L, 1L, 1L), .Dim = c(8L, 2L), .Dimnames = structure(list(         response = c("a", "b", "c", "d", "e", "f", "g", "h"), treatment = c("t1",          "t2")), .Names = c("response", "treatment")), class = "table")


title: Training with very few positives

text: I have a binary classification problem where the fraction of positives is very low, e.g. 20 positives in 10,000 examples (0.2%) **What is an appropriate cross validation scheme for training a classifier with very few positives?** I currently have the following setup:               library(caret)     tmp <- createDataPartition(Y, p = 9/10, times = 3, list = TRUE)     myCtrl <- trainControl(method = "boot", index = tmp, timingSamps = 2, classProbs = TRUE, summaryFunction = twoClassSummary)          RFmodel <- train(X,Y,method='rf',trControl=myCtrl,tuneLength=1, metric="ROC")     SVMmodel <- train(X,Y,method='svmRadial',trControl=myCtrl,tuneLength=3, metric="ROC")     KNNmodel <- train(X,Y,method='knn',trControl=myCtrl,tuneLength=10, metric="ROC")     NNmodel <- train(X,Y,method='nnet',trControl=myCtrl,tuneLength=3, trace = FALSE, metric="ROC")      but I am not getting good performance (my ROC values are < 0.7 for all the classifiers above)


title: Analyzing test scores with different size groups

text: I'm a bit new to statistics; I know some general info about ANOVA, but I'm wondering if it's applicable in the following situation. Let's say I have a group of people who have taken a test, and the group is broken down into different races. However, the number of people of each race is different (30% Caucasian, 50% Asian, 20% African American). I measure the scores on a test for each person, and I want to see if the average score on the test is correlated with race. I know ANOVA is generally used to compare means of different groups, but I'm not entirely sure if it's applicable in this situation when the racial groups are of different size. What kind of analysis do you think I could do to answer this question?


title: Scale parameters -- How do they work, why are they sometimes dropped?

text: I'm having difficulty wrapping my head around scale parameters. How exactly do they work? Why are they sometimes ignored? (in other words, when is it important to preserve them in a calculation?) Thank you for any simple, lay explanations. * * * UPDATE:   What motivated the question was the Gamma distribution; in reading of it with relation to Poisson priors, I saw it described as simply $\propto x^{\alpha-1} e ^{-x /\beta }$ However, the question is more of a general nature. While I get the concept of proportionality, I do not yet fully grasp when it is important for items to be equal versus simply proportional; or rather why so often in the literature I see items referred to as proportional with the scaling factor dropped.


title: What distributions could help describe my uncertainty about a probabilistic forecast?

text: I'm dealing with binary events and I've got people guestimating the chance that they occur. I'd like to translate someone else's guestimate into a probability distribution representing my belief about the true probability. This is called a posterior distribution right? basic bayesian stats? I imagine I want a distribution which looks roughly normal and fairly wide for a forecast of 0.5 but becomes increasingly skewed for forecasts closer to the tails. At the moment I don't have data, I'm just doing exploratory simulations. I'm looking for an R function to that takes a probability forecast and outputs a probability distribution with reasonable parameters. **Added** : For my purposes, even a very very rough approximation will do. Its perfectly fine if statistically literate people roll their eyes. But R code to get me started would be greatly appreciated. **Added** : I fear I still wasn't clear enough about my weird situation. My "expert" is actually only telling me their single best guess for the binary event occurring. I've seen (but not recorded) a bunch of the expert's past best guesses so I have an intuitive sense for what the true distribution of the probability of the events are conditional on an expert's best guess. So I guess I need to play around with beta distribution parameters. Find the parameters that I like for a 0.5 forecast, a 0.25 forecast, 0.1, 0.05 and 0, and then find some smooth functions to interpolate my guesses at what the parameters alpha and beta should be for the space between?


title: Number of features and accuracy/f-measure - Who can explain this results?

text: I'm performing a binary sentiment classification (positive/negative) based on a Naive Bayes classificator and a SVM. To select top k features I use the MRMR algorithm. The model is trained using a 10-fold cross validation and then tested on unseen data. As you can see, I plotted the top k features from 25 to 350 and the accuracy as well as the f-measure since the algorithms have trouble to predict positive classes. Who can help me to understand them? Why do they converge and diverge? And what is the "best" model (main criterion is f-measure) and the best size of k according to the graphs? ![enter image description here](http://i.stack.imgur.com/ZjQ22.png) ![enter image description here](http://i.stack.imgur.com/X1dB0.png) ![enter image description here](http://i.stack.imgur.com/lT1CA.png) ![enter image description here](http://i.stack.imgur.com/bvNpt.png)


title: On-line pearson calculator for students

text: I've been helping a friend with his stats homework. We are now working on Pearson's Correlation, whose calculation is done by a rigorously constructing a table of $x$, $y$, $(x_i-\overline{x})$, $(y_i-\overline{y})$, $(x_i-\overline{x})^2$, $(y_i-\overline{y})^2$, $\sum{(x_i-\overline{x})(y_i-\overline{y})}$ and summing the data to the final formula. There are quite a few Pearson calculators online (e.g. Woflram's), but they all provide a final numeric answer. Is there an online calculator that provides a detailed calculation as well?


title: Time series forecast in R with yearly frequency

text: I have a time series with daily observations over the course of multiple years (interest in topic "superbowl" over time). The seasonality in the data is yearly as well and it is very spiky (almost nothing all year and big increase/spike in January/February). I have started using R for this task (`forecast` package) and have little experience with statistics.               x <- ts(myts, frequency=365)     fit <- HoltWinters(x)     plot(forecast(fit))      This works great and captures the seasonality of the data. Now, I have read more about exponential smoothing (at http://otexts.com/fpp/7/) and understood that the HoltWinters model is one instance of the state space models implemented in ets. Unfortunately, I could not use ets so far since it complains about the high data frequency. I definitely need daily forecast (on the order of 30-60 steps).               fit <- ets(x, 'AAA')     Error in ets(x, "AAA") : Frequency too high      Why can HoltWinters deal with this but not ets? Is there a good workaround? I have the same problem for seasonal ARIMA models and considered splitting up the data in years and using past years as exogenous input. On a side note: How do you usually deal with leap days that screw up your 365 day period? Simply delete them? Thank you very much! PS: I am aware of this: http://robjhyndman.com/researchtips/longseasonality/ However, I couldn't get it too work well on my data, yet. On the other hand, HoltWinters worked fairly well. * * * Thanks for all the helpful comments and discussion. I uploaded the data at http://timalthoff.de/data/data.zip The plot below shows Super_bowl.dat. I took the liberty of including more time series if you'd like to check out more examples. At certain points in time I want to forecast the time series on the order of 60 days. These points in time usually are on the left flank of a big spike that represents a sudden interest in a topic. See example.png for an example (the vertical red lines are these points in time to start an out-of-sample forecast). For more info check out the README. ![enter image description here](http://i.stack.imgur.com/iKy6A.jpg)


title: Standard error of random effects in R (lme4) vs Stata (xtmixed)

text: Please consider this data:               dt.m <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), occasion = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("g1", "g2"), class = "factor"),     g = c(12, 8, 22, 10, 10, 6, 8, 4, 14, 6, 2, 22, 12, 7, 24, 14, 8, 4, 5, 6, 14, 5, 5, 16)), .Names = c("id", "occasion", "g"), row.names = c(NA, -24L), class = "data.frame")      We fit a simple variance components model. In R we have:               require(lme4)     fit.vc <- lmer( g ~ (1|id), data=dt.m )      Then we produce a caterpillar plot:               rr1 <- ranef(fit.vc, postVar = TRUE)     dotplot(rr1, scales = list(x = list(relation = 'free')))[["id"]]      ![Caterpillar plot from R](http://i.stack.imgur.com/fyV1H.jpg) Now we fit the same model in Stata. First write to Stata format from R:               require(foreign)     write.dta(dt.m, "dt.m.dta")      In Stata               use "dt.m.dta"     xtmixed g || id:, reml variance      The output agrees with the R output (neither shown), and we attempt to produce the same caterpillar plot:               predict u_plus_e, residuals     predict u, reffects     gen e = u_plus_e – u     predict u_se, reses          egen tag = tag(id)     sort u     gen u_rank = sum(tag)          serrbar u u_se u_rank if tag==1, scale(1.96) yline(0)      ![enter image description here](http://i.stack.imgur.com/hQKEP.jpg) Clearty Stata is using a different standard error to R. In fact Stata is using 2.13 whereas R is using 1.32. From what I can tell, the 1.32 in R is coming from               > sqrt(attr(ranef(fit.vc, postVar = TRUE)[[1]], "postVar")[1, , ])      [1] 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977      though I can't say I really understand what this is doing. Can someone explain ? And I have no idea where the 2.13 from Stata is coming from, except that, if I change the estimation method to maximum likelihood:               xtmixed g || id:, ml variance      ....then it seems to use 1.32 as the standard error and produce the same results as R.... ![enter image description here](http://i.stack.imgur.com/apbDd.jpg) .... but then the estimate for the random effect variance no longer agrees with R (35.04 vs 31.97). So it seems to have something to do with ML vs REML: If I run REML in both systems, the model output agrees but the standard errors used in the caterpillar plots don't agree, whereas if I run REML in R and ML in Stata, the caterpillar plots agree, but the model estimates do not. Can anyone explain what is going on ?


title: The use of effect size

text: I was able to calculate the sample size using two different effect sizes 0.3 and 0.5 and each gave me different sample sizes. 0.3 gave me 82 and 0.5 gave me 26. I have 45 students in that gradelevel. Does this mean that I am able to use these number of students. What exactly is effect size and why does it affect the sample size. Also, does anyone know how I should write this into my proposal. Thank you very much to Stephan Kolassa and Dandar, their answer to my first question helped me get here.


title: Standard deviation of a dataset comprised of triplicates (industrial stack sampling data)

text: I have 4 triplicate data values. each consist of 3 values, for example, this could be 1:               Concentration:     45     78     66      These data have a mean of `63` and a st.dev. of `16.7` Now I have 4 similar triplicates, from multiple tests, for example:               test:   a     b     c     d             45    56    34    67             78    23    67    34             66    23    56    23          avg     63    34    52.3  41.3      I need the standard deviation of the overall distribution from the 4 triplicates. What is the method to do this? I could average the 4 concentrations, and determine the pooled variance, but would this not be the estimate of the standard deviation if another triplicate test was conducted and the result was the average value? I could determine the standard deviation of the 4 values, but this ignores the fact that these values consist of a triplicate set of data. I need to know what the proper method of determining the standard deviation is.


title: Confidence Interval for a Percent Change

text: Good Morning, A mentor has asked me to calculate the 95% confidence interval surrounding a percent change. To help frame the question, I have calculated the incidence of a disease within the population as follows: **Period 1:** # with disease = 22 population at risk = 2836131 incidence per 100,000 = 0.78 95% confidence interval = 0.49, 1.17 **Period 2:** # with disease = 46 population at risk = 8011777 incidence per 100,000 = 0.57 95% confidence interval = 0.42, 0.77 For reference, I calculated the above incidence rates and their confidence intervals by assuming a Poisson distribution in Stata: cii 2836131 22 , poisson cii 8011777 46 , poisson The end result is that I have reported the decline in the incidence of disease (0.78/100,000 to 0.57/100,000), which is a 26% decrease. My mentor would like to see a 95% confidence interval surrounding the percent change (26%)... If anyone could help by providing a reference or a suggestion of where to go with this, I would be greatly appreciative. Thanks!


title: adjust rating based on number of experiences

text: Suppose I am using a restaurant rating as an explanatory variable in a regression. The rating is defined as $R=\frac{G}{G+B+N+S}$, where $G$ is good, $B$ is bad, $N$ is neutral and $S$ is silent. I have two conceptual issues. First, I want to adjust the rating when it's based on a small number of experiences, perhaps by shrinking it towards the overall mean. Second, people seem overly reluctant to rate highly-attended restaurants ($S$ is very high for places like the local joint that's been around for 30 years). If I plot ratings against $\ln(\text{experiences})$, I get an inverted U shape. Are there any transformations that I can use to remedy these two issue?


title: Principal component analysis, bootstrap and probability of eigenvalue collision?

text: This is really a side project of mine ... while writing on a paper on something totally different! I read (part of ) the excellent paper "FINITE SAMPLE APPROXIMATION RESULTS FOR PRINCIPAL COMPONENT ANALYSIS: A MATRIX PERTURBATION APPROACH" BY BOAZ NADLER This shows among other things that misbehavior of PCA is connected to "crossings of eigenvalues". If you are interested in estimating the eigenvector corresponding to the largest (model) eigenvalue, then if the largest sample eigenvalue frequently comes from a smaller (model) eigenvalue, which is a "crossing", the corresponding sample eigenvector do not correspond to what we try to estimate, and PCA breaks down. So, given a sample, it would be interesting to estimate (for example, using the bootstrap) the probability of crossing. How can we do that? Write the sample covariance matrix as $S$ and a covariance matrix corresponding to a perturbation of the sample $S^*$. (The perturbation could come from the bootstrap, or, conceivably, some other source). The problem now is that the eigensolutions of $S$ and $S^*$, computed by some algorithm, not necessarily are "the same", so we cannot pair of eigenvalues from $S$ with eigenvalues from $S^*$. ¿What does this mean? Write $t \mapsto t S + (1-t) S^* = S_t$. To be able to "pair of" eigenvalues we need an eigensolution of $S_t$ that depends continuously on $t$! The usual algorithms do not give us that. So, my question is: Anybody knows about such algorithms? any other related ideas? papers? Thank you!


title: Predicting the outcome of the next coin toss?

text: Tossing a fair coin is i.i.d. Let's assume that I have that:               Coin = {H,H,H,T,T}      How would you guess the next coin and why? **PS:** I looked for Bernoulli Distribution but couldn't find a logical answer.


title: Sanity check to determine if 'random' data on disk drive may have been tampered with

text: I've got some code that does a Dept of Defense secure erase of disks, and would like to know about a good sanity check to run after the fact. Bottom line, the last write phase of a secure erase of a HDD is to write random data on every byte on the HDD, per the spec. The official standard does not put any constraints on it beyond that. So I use a decent cryptographic algorithm to write data ... but that is not an issue here. What I would really like to determine is what a good sanity check is on a HDD that has been erased. As the spec does not allow one to add any meaning to what is written on the disk, is there a way to assess the probability that the data is random? I read every byte on the surface and create 256 buckets for each possible byte value, and currently count the number of consecutive bytes, and also look at standard deviation to see if any of the byte values are significantly different than others. With that in mind, what is a reasonable # of consecutive bytes for a HDD of X bytes to consider that results are 'inconclusive', or the data might have been tampered with? Same for std deviation, and is there a better model? (Yes, I know somebody could still hide data all over the place if they wanted to, I'm just looking for a way to detect if the process may not have completed, or perhaps the very beginning or end of the disk were never randomized)


title: Are the sample range and sample variance independent when population is normally distributed?

text: If a population is normally distributed, the sample mean and sample variance are independent. What about the sample range and sample variance? Are they independent too? I am trying to derive Studentized Range Distribution, and it might be simpler if it is... or not...


title: calculating variance of sum of predicted costs from two-part model

text: I have a question regarding the calculation of the variance of the sum of predicted health care costs using a two-part regression model. Details are below, but it boils down to how to calculate the covariance for a random variable which is itself the product of two random variables. The two-part regression model consists of a logistic regression component for estimating the probability $p$ of cost > 0 and an OLS component predicting the cost $y$ given $y > 0$. The predicted cost $C_i$ for a given member $i$ is then the product $C_i=p_iy_i$. We can assume that $p_i$ and $y_i$ are independent. I'm calculating the variance of the sum of costs $S=\Sigma_{i=1}^mC_i$ for members $i=1$ to $m$ belonging to a particular provider, assuming member level predicted costs are _not_ independent given they're sharing the same two-part regression model. So running through the math: $Var(S)=\Sigma_{i=1}^m Var(C_i)+\Sigma\Sigma_{i\neq j}^mCov(C_i,C_j)$ where $Var(C_i)$ is: $Var(C_i) = E(C_i^2) - [E(C_i)]^2$ $=E(p_i^2y_i^2)-[E(p_i)E(y_i)]^2=E(p_i^2)E(y_i^2)-[E(p_i)E(y_i)]^2$ $=(Var(p_i)+[E(p_i)]^2)(Var(y_i)+[E(y_i)]^2)-[E(p_i)E(y_i)]^2$ and $Cov(C_i,C_j)$: is $Cov(C_i,C_j)=E(C_iC_j)-E(C_i)E(C_j)$. My first inclination was to use a multivariate Taylor expansion to approximate the term $E(C_iC_j)$. However I run into trouble with a first-order expansion: $Z=g(C_i,C_j)\approx g(\mu_{C_i},\mu_{C_j}) + (C_i-\mu_{C_i})\frac{\delta g(\mu_{C_i},\mu_{C_j})}{\delta C_i} + (C_j-\mu_{C_j})\frac{\delta g(\mu_{C_i},\mu_{C_j})}{\delta C_j}$ since $E(Z)\approx g(\mu_{C_i},\mu_{C_j})=E(C_i)E(C_j)$ which means $Cov(C_i,C_j)\approx 0$, which can't be right. Is there another method for approximating the covariance of $C_i$ and $C_j$? Or did I make a mistake in my calculations? Thanks!


title: Transform non-normal data to normality by rescoring columns

text: I have a vector with fluency of translation - (0, 0, 1, 3, 3, 3 ,3 ....) The problem is that it is made by people (for example someone gives too much 3 but only a bit of 2) and we want to normalize them. We have a vector from each person separately. It is possible to do that by only rescoring columns of histogram like a picture? ![fluency](http://i.stack.imgur.com/BNvem.png) There is no matter if the new x-values would be other than from (-3,3). Is there a method how to do that? Or a function in R?


title: How to calculate probability contingent on multiple events?

text: This isn't a homework problem, but something I just thought of when watching a crime show. If I am trying to determine the probability of [insert Conclusion] or C, and I know C is comprised of A being true (which is comprised of A1 and A2 being true, each which have a 50% chance of happening) and B being true (which has a 20% chance of happening), how would I formulate this?


title: One-sample KS test's D-statistic peculiarity in R

text: First of all, let me explain what I'm after. I have a bunch of archaeological sites (point data) in GIS and a DEM-derived slope map. I classified the slope map into 1-degree interval classes and got the values for the sites. Now I wanted to test whether the values of the sites are significantly different than would be their random placement on the slope map, meaning that the slope of the terrain was an important factor when choosing a place for settlement. This seemed to be a good case for 1-sample KS test as the data is continuous and the reference distribution is actually known (the slope map itself). Also, the distributions are clearly non-normal. So, I used `R`'s `ks.test` function so that the numeric vector x contained the classified values for the sites and y was set the cumulative distribution of the background distribution values via `ecdf` function: the classification into 1-degree bands was, of course, the same. However, the result is this:               One-sample Kolmogorov-Smirnov test          data:  slope_sites     D = 0.4046, p-value < 2.2e-16     alternative hypothesis: two-sided      Not only is the p-value ridiculously small, the D statistic is also peculiar. When plotting the two empirical cumulative distribution functions, the difference is much smaller. Not wanting to believe such result I computed the D statistic myself:               max(abs(ecdf(slope_background)(0:55) - ecdf(slope_sites)(0:55)))      The 55 is set as limit as 55 degrees is the maximum recorded value Anyway, the resulting D value was 0.1650732, which seems more realistic and is in accordance to the plotted ECDFs. Any thoughts on that? I suspect I'm doing something fundamentally wrong though. Also, is there any other, perhaps better way to test this kind of hypothesis? I tried `ad.test` with the same parameters and I suspect that p-value of 4.38e-06 is also not exactly what would be expected. The histogram of the site distribution: ![Hist site](http://i.stack.imgur.com/vPAWw.png) Cumulative distribution functions of the sites (red) and background data (black dashed line): ![ECDF](http://i.stack.imgur.com/j5fxk.png)


title: Are the relations in fixed, random and mixed effect models and multilevel models causal?

text: In fixed, random and mixed effect models, and multilevel models, the response random variable is represented as a function of some explanatory variables and random errors. I was wondering if the relations implied by them are considered causal, and therefore used in causal inference? Thanks!


title: R: GLMM using lme4

text: I've got a conditional logistic regression setup using `clogit` in `R` like this:               m<-clogit(PHENO==2 ~ x + as.factor(COVAR[,1]) + strata(COVAR[,2]) )      I wanted to try doing GLMM analysis in `R`. I'm a little confused on the syntax for the `lme4` package in `R`. `COVAR[,2]` is the matching variable in my data. Can someone explain to me the difference between some of these statements:               m <-  lmer(PHENO==2 ~ x + as.factor(COVAR[,1]) + (x|COVAR[,2]) )     m <-  lmer(PHENO==2 ~ x + as.factor(COVAR[,1]) + (1|COVAR[,2]) )      Which one is more appropriate?


title: Comparison of symptom profiles: binomial GLM? (SPSS, R)

text: I have nine depression symptoms in a sample of $n=2000$; symptoms are dichotomized (0= not present, 1= present) and can thus be used as counts. I want to test whether the symptoms differ from each other in their counts, but I want to do it in one multivariate test, because with 9 symptoms I'd otherwise have 36 comparisons to compute. My problem is: against what data do I compare the symptom profile? Do I understand correctly that I would need to simulate data now in which symptoms are equally distributed (using the mean count of the nine symptoms in my data)? To simplify, assume we have the counts $S_1=100$, $S_2=200$, and $S_3=300$ present in my sample. I would now test this against a simulated sample in which $S_1=S_2=S_3=200$? What program would you recommend for this procedure in case it is the one you would favor, and what do I have to be careful about when simulating the counts? Or would you recommend another procedure altogether? Thank you!


title: Explanatory variables with many zeros

text: I am trying to fit a linear model to a price response variable. Many of the predictor variables consist of mainly zeros. For example, one possible predictor variable is "drill holes". Not many parts have a drilled hole, but if they do it would make sense that it affects the price. I am using the `caret` package in `R` to train the model and choose the appropriate variables. I have already removed all variables with zero variance. I have found a lot of literature about count data for a response variable with many zeros and zero-inflation models. But what I am wondering is, how should explanatory variables with many zeros (many are NOT count data) be handled? Is there an appropriate transformation? Or are explanatory variables with many zeros allowable since I am dealing with explanatory variables and not the response variable?


title: Treatment Effects: Manova

text: Suppose we are measuring blood pressure three times a day for two treatments $A$ and $B$. 10 people are getting treatment A and 10 people are getting treatment B. Also suppose we find that there is a treatment and time interaction. To test for equal treatment effects, how would we incorporate the interaction with time? Would we do the following:               proc glm;       class trt;       model bp1 bp2 bp3 =trt;       manova h=trt m=bp1+bp2+bp3;      or would we do the following:                proc glm;           class trt;           model bp1 bp2 bp3 =trt;           manova h=trt m=(bp1+bp2+bp3)/3;      where we divide by the three time points?


title: Accounting for pretreatment differences (linear mixed effects model) using R/nlme

text: I have an experiment that is as follows:   * 4 years of data (1 year pre-treatment, 3 years during treatment)   * 20 plots total   * 5 plots of each type (Control, Treatment A, Treatment B, Treatment A+B)   * Each plots is divided into 4 quadrants in which measurements are recorded. Initial stats revealed pretreatment differences between plot types. The data I would like to use in a linear mixed-effects model include:   1. my response variable;   2. BA - before/after (as "1" or "2");   3. Treatment A (as "yes" or "no");   4. Treatment B (as "yes" or "no");   5. The interaction of treatments and Year;    6. Plot and quadrant as nested random effects. I've created a linear mixed effects model using the nlme package in R that currently looks like this:               model<-lme(response variable ~ BA * treatmentA * treatmentB * year ,                method= "REML", data=data, random=~1|plot/quad, na.action=na.omit)      I get the following error message: > Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at > level 0, block 1 I can run the model successfully without BA or without year but am unable to get results when I run it with both of them. Any advice or suggestions on how to account for pretreatment differences and get my LME to work happily?


title: Why does the p-value of a composite null hypothesis have a supremum attached it?

text: I noticed that there is a definition of the p-value in my textbook. It is defined as the p-value of a composite null hypothesis and it says the following: ![enter image description here](http://i.stack.imgur.com/7Y4yL.jpg) I have no idea why it is written with a supremum. I've spent hours pondering this, does anyone have enough of a background to help me with this? Thank you!


title: Explanatory variables in a Lin-log model

text: I have a data set and I want to fit a `Lin-log` model. Is it possible to apply the log transformation only to some of the explanatory variables or should the log transformation be applied to all of the explanatory variables?


title: Which distribution fits data better?

text: I use `fitdistr` in R to select which distribution fits my data best. I tried `cauchy`, `weibull`, `normal`, and `gamma`. The log-likelihood is $-329.8492$ for `cauchy`, $-277.4931$ for `gamma`, $-327.7622$ for `normal`, $-279.0352$ for `weibull`. Which one is the best? The one with the largest value (i.e., gamma) or the one with the largest absolute value (i.e., Cauchy)?


title: Plotting top 5 most frequent factors using R

text: In R I have a data frame f with headers a, b, c. Column c has factor data. plot(f$c) does this ![the bad plot](http://i.stack.imgur.com/eZyQh.png) How can I limit the plot to the top 5 most frequently occurring factors of column c in descending order? The following works just fine:               library(ggplot2)     data(diamonds) # clarity is a good categorical variable     with(diamonds, barplot(rev(sort(table(clarity))[1:5])))      Is there a neat way to post this on SO after the fact? ![enter image description here](http://i.stack.imgur.com/unieb.png)


title: How to analyse this simulation

text: So, here's the situation: I have some simulation results, from a simulation driven by two parameters. One of these parameters has two settings, the other has four each was replicated ten times for a total of 80 sets of results. Each simulation produces results for each of four groups. The progress of each group was measured by seven variables. Thus there are 28 results from each simulation. I am certain (due to the implementation of the simulation) that these results are to some extent influencing each other. Three of the outcome variables describe how well the group did in intergroup competitions (which influenced the quantity of rescources available to the group) and four describe the behaviours that the agents learned to use (which is influenced by the quantity of resources available to the group) What I'd like to test is: Did the two parameters influence the four behaviour outcomes and if so in what direction? Did the groups that performed well or poorly in the intergroup competitions (as measured by three variables) behave differently (as measured by these four variables) and if so how? For some reason I'm struggling to get my head around the right approach to take here, but I've got the nagging sensation that it should be obvious.


title: X-mean algorithm BIC calculation question

text: I'm having trouble understanding some of the formulas in this paper related to BIC calculation (Dan Pelleg and Andrew Moore, _X-means: Extending K-means with Efficient Estimation of the Number of Clusters_ ). First the variance equation:   * R - number of points   * K - number of clusters   * $\mu_i$ - centroid associated with ith point.   * $\sigma^2 = \frac{1}{R-K}\sum_{i}(x_i - \mu_{(i)})^2 $ The log likelyhood then uses this sigma. Am I reading this right, they're using 1 covariance matrix for all clusters (see quote below, they are)? This makes no sense. If you have 5 clusters, each one is a Gaussian according to k-means algorithm. So wouldn't it make sense to compute covariance $\sigma^2_i$ for each cluster and use that? My second question is regarding number of parameters to use in the BIC score. The paper mentions > Number of free parameters $p_j$ is simply the sum of K-1 class > probabilities, M*K centroid coordinates, and one variance estimate. How do you get the K-1 class probabilities? I could do # of points in class i / total number of points. But then it's K-1, which probability is left out of the sum? P.S. If anyone has a nicer paper on estimating k using similar methods I'd like to read that as well. At this point I'm not too concerned with speed. Thanks for your help.


title: How to model the effect of "words occurring over time"?

text: As time goes on, some new words may occur. But standard Latent Dirichlet Allocation assume a fixed vocabulary. Is there some variation of the Latent Dirichlet Allocation model which can represent this phenomenon? Thanks!


title: Is the probability of a unique DNA match the same, irrespective of the size of the dataset?

text: One can argue that the evidential value of a database hit increases when the database gets larger. Indeed, having a unique hit in a database simply means that all the other persons in the database are innocent. Therefore, if the database is larger we have the extra information that many more people are innocent. Clearly, this increases the weight of the evidence against the suspect.   * Is this reasoning correct?    * **Is the probability of a unique match the same, irrespective of the size of the dataset?** This is how I see it: if the database contains everybody, then a unique hits tells us with certainty who the donor of the DNA prole is. So it really is true that a larger database is stronger evidence. But then I started thinking: suppose the database consists of 16,000,000 proles, and suppose that the DNA prole we are looking for has frequency 1/1,000,000 (one in a million). Is it still reasonable to believe that the underlying model is correct?


title: How do I ensure my implementation of a MLP from Weka is correct?

text: I have been training MLPs for a binary classification problem in the Weka explorer and now have one with an acceptable level of accuracy. I've written some code to parse the text of the model which is output, and write it out as array literals (in my C# code) for the weights and bias terms. I then have my own code to implement the neural network output. _The problem is that my confusion matrix differs from Weka's; for example one of the values was 86.8% in Weka and 91.3% in my code. How can I explain this difference?_ I've downloaded someone else's C# neural network code and their matrix was the same as mine. So I doubt it is because my implementation of the neural network output is incorrect. I've manually checked that the weights and the normalisation ranges are the same. Weka is running 10-fold cross validation, whereas my code is evaluating the error on the whole data set. _Is there a way in Weka to run the classifier against a data set, and obtain the confusion matrix, but without changing the classifier? Or, to run the classifier and get the actual output values_ , so I can compare them with those in my code?


title: How much do answers A,B,C link with the answer on question D?

text: I have 9 questions in 3 different blocks where you answer with a number (1-7). And there is a question (how did you like this place overall?) that is answered (1-10). I need a way to find out how much each block (and maybe also each question) affects the tenth question. Preferably using excel or R. The number of respondents are in the thousands if that makes any difference.


title: How do I do group wise clustering in R?

text: How do I do group wise clustering in R? Hi all, I have N x K data matrix, where N is the number of observations, K is the number of variables. The N observations fall into M categories or groups. Now I want to cluster the groups, instead of the observations, how do I do that? i.e. the clustering would be at the group level... Thanks a lot for your help!


title: Testing the quality simple linear regression model

text: I have made a simple linear regression model based on two sets of monthly data from 1960 to 2008, using only data from 1960 to 2000. I was now wondering how I would go about back-testing this model so I can comment on its quality. So far I have predicted the values of the dependent variable using the independent variable from 2000 to 2008 and made a residual plot. What else can I do to test the quality of my regression model?


title: Compositional data model and aggregated multinomial logit model

text: What is the difference between a compositional data model using the additive log-ratio (ALR) transformation and an aggregated multinomial logit model?


title: How to report two-way anova results?

text: After finding a significant treatment in two-way anova, how do I report where the differences are? Every text I have read has left me at: Ah, the results are significant, so let's move on. In one-way anova I could use Tukey's HSD to find the means that differed, making my report much more effective. Right now I don't know how to do that with two-factor anova. Here is an example from a text done in R: Is there a difference in Friday tardy rates at different plants?               > mlate             day  plant absences     1   march-4 plant1       19     2  march-11 plant1       22     3  march-18 plant1       20     4   march-4 plant2       18     5  march-11 plant2       20     6  march-18 plant2       16     7   march-4 plant3       27     8  march-11 plant3       32     9  march-18 plant3       28     10  march-4 plant4       22     11 march-11 plant4       27     12 march-18 plant4       26          > anova(lm(absences ~ plant+day, data=mlate))     Analysis of Variance Table          Response: absences               Df  Sum Sq Mean Sq F value    Pr(>F)         plant      3 216.250  72.083  41.191 0.0002134 ***     day        2  30.167  15.083   8.619 0.0172128 *       Residuals  6  10.500   1.750                           ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       So I can see that both the plant and the day have an effect on worker's being late, but how do I compare the means like I would with a tukey's hsd in a one- way test?


title: How to deal with multicollinearity without dropping variables

text: In establishing a linear model, if two variables are highly correlated, then it turns out to have high VIFs. If we applied data reduction techniques, like PCA, we will lose explainability to some extent. In my case, we can't find summary scores that can be put into the model instead of individual variables. But for some other reasons, both variables have to be kept in the model. In this case, what should we do to fix this problem?


title: How can I get the CI of a couple of numbers?

text: Is it possible to get the Confidence Interval of a couple of numbers (or more) without knowing the distribution or anything like that ? Thanks.


title: Is this prior distribution exchangeable?

text: I am reading the book Bayesian data analysis written by Gelman. And I do not confirm the answer of this question:   suppose it is known a priori that the $2J$ parameters $\theta_1,\ldots,\theta_{2J}$ are clustered into two groups, with exactly half being drawn from a $N(1,1)$ distribution, and the other half being drawn from a $N(-1,1)$ distribution, but we have not observed which parameters come from which distribution.   $(a)$ Are $\theta_1,\ldots,\theta_{2J}$ exchangeable under this prior distribution?   $(b)$ Show that this distribution cannot be written as a mixture of $iid$ components.   $(c)$ Why can we not simply take the limit as $J\rightarrow \infty$ and get a couterexample to $de \space Finetti's$ theorem? IMO, we have no any information about ordering or grouping of the parameters can be made, so the parameters $\theta_1,\ldots,\theta_{2J}$ are exchangeable in this prior distribution?   Besides, not all parameters are identity distributed, can we simply say this distribution cannot be written as a mixture of $iid$ components?   We cannot take the limit as $J\rightarrow \infty$ and get a couterexample to $de \space Finetti's$ theorem since this distribution is not exchangeable. Right? If you have any idea please write your comment. Let me know I do a right answer or not. Thanks :)


title: Arellano-Bond estimator in Stata

text: I have a question about Arellano-Bond model in Stata (xtabond/xtabond2). The slopes I get, are they for levels or differences of values? My model to be estimated has a form of (D is first difference): DY=a+DX1+DX2+.... So should I use already differentiated variables in my analysis or levels variables to get this? I know A-B first-differiantes everything at first, so if I use already differentiated data, are the slopes for differences of differences? In most examples levels variables are used. Risto from Finland


title: Difference between using a dummy variable in a multiple regression and clustering on that dummy

text: In a multiple regression, the coefficients become insignificant using `vce(cluster year)` in Stata. Will including year dummies solve the problem? How are these two approaches different?


title: Function of an efficient estimator

text: Let T be an efficient estimator for a parameter $\theta$. Let $g(\cdot)$ be a (real valued) non-linear function. Is $g(T)$ an UMVUE for $g(\theta)$?


title: What is the connection between credible regions and Bayesian hypothesis tests?

text: In frequentist statistics, there is a close connection between confidence intervals and tests. Using inference about $\mu$ in the $\rm N(\mu,\sigma^2)$ distribution as an example, the $1-\alpha$ confidence interval $$\bar{x}\pm t_{\alpha/2}(n-1)\cdot s/\sqrt{n}$$ contains all values of $\mu$ that aren't rejected by the $t$-test at the significance level $\alpha$. Frequentist confidence intervals are in this sense inverted tests. (Incidentally, this means that we can interpret the $p$-value as the smallest value of $\alpha$ for which the null value of the parameter would be included in the $1-\alpha$ confidence interval. I find that this can be a useful way to explain what $p$-values really are to people who know a bit of statistics.) Reading about the decision-theoretic foundation of Bayesian credible regions, I started to wonder whether there is a similar connection/equivalence between credible regions and Bayesian tests. Is there a general connection? If there is no general connection, are there examples where there is a connection? If there is no general connection, how can we see this?


title: What is the quantity "risk" in survival analysis?

text: From Wikipedia > ... the hazard ratio is a relative measure of effect and tells us nothing > about **absolute risk**. What is the definition of "absolute risk" in survival analysis? Does it depend on time? Is it the value of the survival function at some time? Thanks!


title: Recommendations for first time teacher (Intro to Biostatistics)

text: I'm teaching my first class this fall (Intro to Biostatistics). Anyone have any suggestions for teaching statistics better? Perhaps some example you wish your first teacher had used? I'm using Principles of Biostatistics by Pagano and Gauvreau. **EDIT: DETAILS** This class is an online class that meets twice a week for 1.5 hours. Students will listen to me lecture while viewing a powerpoint/beamer presentation (boring?) mixed with a little tablet/pen action (exciting?) This class is VERY beginner statistics stuff, taught to mostly biomedical engineers (undergrads) and some non-statistics graduate students (nursing, medical students, public health, etc) **Syllabus** : 1) What is biostatistics? 2) Probability 3) Diagnostic tests (i.e. specificity, sensitivity, ROC curves. Mostly here b/c it allows us to apply some of the things we learned about in probability e.g. Bayes rule) 4) Distributions 5) Sampling Distributions 6) Confidence Intervals 7) Hypothesis testing (one sample, two sample, proportions) 8) Power and Sample size calculations 9) Non-parametric methods 10) Contigency tables (chi-square test, fishers test, mcnemars test, relative risk, odds ratios 11) correlation Main object is for students to learn the core concepts of statistical inference e.g how do you quantify the question "Which drug is better?" stuff like that. For the sections listed above, do you have any words of advice/caution for teaching them. **For example** : I've seen/heard of many different approaches to introducing students to the wald test vs the score test when doing inference on proportions. If taught poorly, students easily get confused ("Why are there two of these?", "Which one do I use?", "They look the same to me.") Some teachers don't even mention these names and just say: Do this for confidence intervals and this other thing for hypothesis testing. How would you approaches this problems or others like it?


title: How to think of features in NLP problems

text: I am working on a Named Entity Recognition (NER) project. Instead of using an existing library, I decided to implement one from scratch because I wanna learn the basics of how PGMs work under the hood. I converted the words in sentences into feature vectors. The features are manually picked by me, and I can only think of roughly ~20 features (such as: "Is the token capitalized?", "Is the token an English word?", etc.). However, I've heard good NER algorithms represent tokens using way more than 20 features, sometimes hundreds of features. How do they manage to think of so many features? Are there any recommended best practices in feature construction?


title: Comparing observed data to simulated data?

text: I have a set of observed values (shown with black dots in the figure) that I would like to compare to some simulated data (100 simulated datasets shown as box plots with quartiles, extremes (excluding outliers) as whiskers, and outliers as white dots). ![Box plots](http://i.stack.imgur.com/UARxG.png) The observed values are outside of the 95% confidence intervals of the simulated dataset so it seems pretty obvious that there is a significant difference between the observed and the simulated; but in some cases, the difference might not be so obvious. Which statistical test (preferably in R) can I do on this data to get a P-value for the significant difference between the observed and the simulated sets?


title: How can I determine accuracy of past probability calculations?

text: I do not study statistics but engineering, but this is a statistics question, and I hope you can lead me to what I need to learn to solve this problem. I have this situation where I calculate probabilities of 1000's of things happening in like 30 days. If in 30 days I see what actually happened, how can I test to see how accurately I predicted? These calculations result in probabilities and in actual values (ft). What is the method for doing this? Thanks, CP


title: SEM: Do I still have to look at the fit of measurement models when overall fit is good?

text: I have a question concerning structural equation modeling: If the overall fit of my model is very good like in the following model, do I still have to test whether the fit of each measurement model is good as well? And if yes, why is that? And what can I do if my whole model includes measuresment models, that just cannot be tested separately because there are only 2 indicators for each latent variable? Thank you in advance for your help! ![SEM](http://i.stack.imgur.com/XhGJt.jpg)


title: Confidence intervals for empirical CDF

text: I have 100 data points from a random process. How would I go about placing a confidence interval around the estimate of $\Pr(X>x)$? The distribution function is unknown and positively skewed. My first inclination would be to use a bootstrap based on the material I have read for this class, but is there some other way to do this?


title: exact Monte Carlo simulation of a large sample histogram

text: (Asked this on MathOverflow, got redirected here) Say I want to create a histogram of $N$ samples from some simple compactly supported distribution on $\mathbb{R}$, where $N$ is very large, say $N = 10^{30}$. The histogram has $K$ disjoint bins, where $K$ is a more reasonable number like $10$ or $1000$. Obviously it's not feasible for me to directly draw $N$ Monte Carlo samples from my distribution and bin them up to form the histogram. However, it seems to me that there might be some correct method which works by sampling the number of counts in each bin, one at a time, for a total of only $K$ samples. I'm looking for some help finding such a method. The number of counts in each bin is a binomial random variable with parameters that can easily be calculated by integrating the distribution over the bin interval. So if I have a good way to simulate binomial random variables with large means, I can simulate the number of counts in each bin using only $K$ Monte Carlo draws of a binomial random variable. The problem is that the total counts in the bins are correlated by the constraint that they must add up to $N$. My method will produce a random number of total counts which will almost certainly not be $N$. I can think of a couple more sophisticated methods that would avoid this problem, but they create thornier ones - and the bottom line is, I don't know how to prove that any of these heuristic methods are correct. Can anybody think of an $O(K)$ algorithm that generates a provably correct (or provably nearly correct) histogram sample for this kind of problem? More formally, a correct method for sampling the random vector $H \in \mathbb{Z}^K$ whose entries are the histogram counts? If not, what's the best that I can hope for?


title: Multicollinearity among spatial predictors

text: I am trying to model the effect that proximity to transportation infrastructure has on home prices. It is common (and reasonable) to use some type of spatial dependence model, such as the Spatial Durbin Model, $$y = \rho W y + X\beta + XW\lambda + \varepsilon $$ In this model, the outcome is a result of   * autocorrelation $\rho W y$ (neighbor's home prices)    * individual effects $X\beta$ (own home characteristics)    * and neighbor effects $XW\lambda$ (neighbor's characteristics) My question: in the case of a spatially-dependent explanatory variable (such as how far each observation is from a transit station), does this model introduce collinearity between the $X\beta$ and $XW\lambda$ terms? After all, if I live two miles away from the nearest train station, my neighbors by definition also live about two miles from the same train station. If there is multicollinearity, how will this effect the $\beta$ and $\lambda$ estimates? Will they be consistent/unbiased? If so, how might I correct for this?


title: Parametric, semiparametric and nonparametric bootstrapping for mixed models

text: Following grafts are taken from this article . I'm newbie to bootstrap and trying to implement the parametric, semiparametric and nonparametric bootstrapping bootstrapping for linear mixed model with `R boot` package. ![enter image description here](http://i.stack.imgur.com/b24bi.png) ![enter image description here](http://i.stack.imgur.com/2fTO0.png) ![enter image description here](http://i.stack.imgur.com/zANep.png) ![enter image description here](http://i.stack.imgur.com/VfGOQ.png) ![enter image description here](http://i.stack.imgur.com/Y2Ng0.png) **R Code** Here is my `R` code:               library(SASmixed)     library(lme4)     library(boot)          fm1Cult <- lmer(drywt ~ Inoc + Cult + (1|Block) + (1|Cult), data=Cultivation)     fixef(fm1Cult)               boot.fn <- function(data, indices){      data <- data[indices, ]      mod <- lmer(drywt ~ Inoc + Cult + (1|Block) + (1|Cult), data=data)      fixef(mod)      }          set.seed(12345)     Out <- boot(data=Cultivation, statistic=boot.fn, R=99)     Out      **Questions**   1. How to do parametric, semiparametric and nonparametric bootstrapping for mixed models with `boot` package?   2. I guess I'm doing nonparametric bootstrapping for mixed model in my code. I found these slides but could not get the R package `merBoot`. Any idea where I can get this package. Any help will be highly appreciated. Thanks in advance for your help and time.


title: What test can I use to prove there is an upward trend in time series of ratios assesment in R?

text: I have a database where experiments are stored. Each experiment has a flag, which says that it was properly registered. I am trying to show with a statistical test, that there is a trend in time that experients are increasingly properly registered (ratio column; rising from 56% in 2006 to 72% in 2011)               > print(ascii(c,include.rownames=F,format='nice'),type='org')     | year | all_experiments | properly registered | ratio |     |------+-----------------+---------------------+-------|     | 2006 | 6431            | 3604                | 56.04 |     | 2007 | 7013            | 3990                | 56.89 |     | 2008 | 8285            | 4899                | 59.13 |     | 2009 | 7523            | 5063                | 67.3  |     | 2010 | 7296            | 5210                | 71.41 |     | 2011 | 7243            | 5243                | 72.39 |      What statistical approach should I use to prove it.   1. Do I do take just the percentages, plot them and show that I can fit a line through the data and point to the slope of that line: how do I do that? Is it enough?    2. Is there any more complicated test which takes into into account the numbers behind the percentages, and how would I do it? I also need to compute the analysis within R so a generic help is useful but also hints to what commands to use within R would be best. The data partitioned by year come from the same dataset and were produced like this:               a = ddply(s1b,.(format(First.Received, "%Y")) ,.fun=function(x) nrow(x))     b = ddply(s1b,.(format(First.Received, "%Y")) ,.fun=function(x) nrow(subset(x,properly_registered==1)))     ...


title: How do I figure out what kind of distribution represents this data on ping response times?

text: i've sampled a real world process, network ping times. The "round-trip-time" is measured in milliseconds. Results are plotted in a histogram: ![alt text](http://i36.tinypic.com/2qtfp6f.jpg) Ping times have a minimum value, but a long upper tail. i want to know what statistical distribution this is, and how to estimate its parameters. Even though the distribution is not a normal distribution, i can still show what i am trying to achieve. The normal distribution uses the function: ![alt text](http://upload.wikimedia.org/math/9/e/1/9e1e4a3af93c9680ba75669a0b69fbf6.png) with the two parameters   * μ (mean)   * σ2  (variance) ## Parameter estimation The formulas for estimating the two parameters are: ![alt text](http://upload.wikimedia.org/math/3/2/4/324081bc08a66fbeba64f43a3dbd5e1a.png) Applying these formulas against the data i have in Excel, i get:   * μ = 10.9558 (mean)   * σ2  = 67.4578 (variance) With these parameters i can plot the " _normal_ " distribution over top my sampled data: ![alt text](http://i34.tinypic.com/2mqu3wl.jpg) Obviously it's not a normal distribution. A normal distribution has an infinite top and bottom tail, and is symmetrical. This distribution is not symmetrical. * * * >   * What principles would i apply, what flowchart, would i apply to > determine what kind of distribution this is? >   * Given that the distribution has no negative tail, and long positive > tail: what distributions match that? >   * Is there a reference that matches distributions to the observations > you're taking? > And cutting to the chase, what is the formula for this distribution, and what are the formulas to estimate its parameters? * * * i want to get the distribution so i can get the "average" value, as well as the "spread": ![alt text](http://i38.tinypic.com/2i27sp4.jpg) i am actually plotting the histrogram in software, and i want to overlay the theoretical distribution: ![alt text](http://i38.tinypic.com/nqpnck.jpg) **Note:** Cross-posted from math.stackexchange.com * * * **Update** : 160,000 samples: ![enter image description here](http://i.stack.imgur.com/k51W5.png) Months and months, and countless sampling sessions, all give the same distribution. There _must_ be a mathematical representation. * * * Harvey suggested putting the data on a log scale. Here's the probability density on a log scale: ![enter image description here](http://i.stack.imgur.com/2gHw9.png) **Tags** : sampling, statistics, parameter-estimation, normal-distribution * * * It's not an answer, but an addendum to the question. Here's the distribution buckets. i think the more adventurous person might like to paste them into excel (or whatever program you know) and can discover the distribution. The values are normalized               Time    Value     53.5    1.86885613545469E-5     54.5    0.00396197500716395     55.5    0.0299702228922418     56.5    0.0506460012708222     57.5    0.0625879919763777     58.5    0.069683415770654     59.5    0.0729476844872482     60.5    0.0508017392821101     61.5    0.032667605247748     62.5    0.025080049337802     63.5    0.0224138145845533     64.5    0.019703973188144     65.5    0.0183895443728742     66.5    0.0172059354870862     67.5    0.0162839664602619     68.5    0.0151688822994406     69.5    0.0142780608748739     70.5    0.0136924859524314     71.5    0.0132751080821798     72.5    0.0121849420031646     73.5    0.0119419907055555     74.5    0.0117114984488494     75.5    0.0105528076448675     76.5    0.0104219877153857     77.5    0.00964952717939773     78.5    0.00879608287754009     79.5    0.00836624596638551     80.5    0.00813575370967943     81.5    0.00760001495084908     82.5    0.00766853967581576     83.5    0.00722624372375815     84.5    0.00692099722163388     85.5    0.00679017729215205     86.5    0.00672788208763689     87.5    0.00667804592402477     88.5    0.00670919352628235     89.5    0.00683378393531266     90.5    0.00612361860383988     91.5    0.00630427469693383     92.5    0.00621706141061261     93.5    0.00596788059255199     94.5    0.00573115881539439     95.5    0.0052950923837883     96.5    0.00490886211579433     97.5    0.00505214108617919     98.5    0.0045413204091549     99.5    0.00467214033863673     100.5   0.00439181191831853     101.5   0.00439804143877004     102.5   0.00432951671380337     103.5   0.00419869678432154     104.5   0.00410525397754881     105.5   0.00440427095922156     106.5   0.00439804143877004     107.5   0.00408656541619426     108.5   0.0040616473343882     109.5   0.00389345028219728     110.5   0.00392459788445485     111.5   0.0038249255572306     112.5   0.00405541781393668     113.5   0.00393705692535789     114.5   0.00391213884355182     115.5   0.00401804069122759     116.5   0.0039432864458094     117.5   0.00365672850503968     118.5   0.00381869603677909     119.5   0.00365672850503968     120.5   0.00340131816652754     121.5   0.00328918679840026     122.5   0.00317082590982146     123.5   0.00344492480968815     124.5   0.00315213734846692     125.5   0.00324558015523965     126.5   0.00277213660092446     127.5   0.00298394029627599     128.5   0.00315213734846692     129.5   0.0030649240621457     130.5   0.00299639933717902     131.5   0.00308984214395176     132.5   0.00300885837808206     133.5   0.00301508789853357     134.5   0.00287803844860023     135.5   0.00277836612137598     136.5   0.00287803844860023     137.5   0.00265377571234566     138.5   0.00267246427370021     139.5   0.0027472185191184     140.5   0.0029465631735669     141.5   0.00247311961925171     142.5   0.00259148050783051     143.5   0.00258525098737899     144.5   0.00259148050783051     145.5   0.0023485292102214     146.5   0.00253541482376687     147.5   0.00226131592390018     148.5   0.00239213585338201     149.5   0.00250426722150929     150.5   0.0026288576305396     151.5   0.00248557866015474     152.5   0.00267869379415173     153.5   0.00247311961925171     154.5   0.00232984064886685     155.5   0.00243574249654262     156.5   0.00242328345563958     157.5   0.00231738160796382     158.5   0.00256656242602444     159.5   0.00221770928073957     160.5   0.00241705393518807     161.5   0.00228000448525473     162.5   0.00236098825112443     163.5   0.00216787311712744     164.5   0.00197475798313046     165.5   0.00203705318764562     166.5   0.00209311887170926     167.5   0.00193115133996985     168.5   0.00177541332868196     169.5   0.00165705244010316     170.5   0.00160098675603952     171.5   0.00154492107197588     172.5   0.0011150841608213     173.5   0.00115869080398191     174.5   0.00107770703811221     175.5   0.000946887108630378     176.5   0.000853444301857643     177.5   0.000822296699600065     178.5   0.00072885389282733     179.5   0.000753771974633393     180.5   0.000766231015536424     181.5   0.000566886361087923


title: Gaussian process predictor

text: I am building GP regressor , my input data is 1-d column vector and so is my target. I have divided my data into training and testing sets. I trained the model to learn the hyper-paramters and then plugged in my input test cases, and I expected to see their corresponding outputs as in the following lines of code :               % Split the data     XTrain = BP(1:193,1);     YTrain = BP(1:193,2);     XTest =  BP(194:end,1);     YTest =  BP(194:end,2);     % set meand, cov , lik functions and hyp     meanfunc = {@meanSum, {@meanLinear, @meanConst}};     covfunc = {@covMaterniso, 3};     likfunc = @likGauss;     hyp.mean = [0.5; 1];     hyp.cov = log([1/4; 1]);      hyp.lik = log(0.1);     %Training and learning parameters      hyp = minimize(hyp, @gp, -100, @infExact, meanfunc, covfunc, likfunc, XTrain, YTrain);     % Compute predictions      [m s fm fs] = gp(hyp, @infExact, meanfunc, covfunc, likfunc, XTrain, YTrain, XTest);     % Code ends here       Is what I am doing right? and how can I interpret the mean?


title: Are degrees of freedom in lmerTest::anova correct?

text: I am analyzing the results of a reaction time experiment in R. I ran a repeated measures ANOVA (1 within-subject factor with 2 levels and 1 between-subject factor with 2 levels). I ran a similar linear mixed model and I wanted to summarize the lmer results in the form of ANOVA table using `lmerTest::anova`. Don't get me wrong: I did not expect the identical results, however I am not sure about the degrees of freedom in `lmerTest::anova` results. ~~It seems to me it rather reflects an ANOVA with no aggregation on subject-level.~~ I am aware of the fact that calculating degrees of freedom in mixed-effect models is tricky, but `lmerTest::anova` is mentioned as one possible solution in the updated `?pvalues` topic (`lme4` package). Is this calculation correct? Do the results of `lmerTest::anova` correctly reflect the specified model? **Update:** I made the individual differences larger. The degrees of freedom in `lmerTest::anova` are more different from simple anova, but I am still not sure, why they are so large for the within-subject factor/interaction.               # mini example with ANT dataset from ez package     library(ez); library(lme4); library(lmerTest)          # repeated measures ANOVA with ez package     data(ANT)     ANT.2 <- subset(ANT, !error)     # update: make individual differences larger     baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)     ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]          anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum),        within = .(direction), between = .(group))     anova.ez          # similarly with lmer and lmerTest::anova     model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)     lmerTest::anova(model)          # simple ANOVA on all available data     m <- lm(rt ~ group * direction, data = ANT.2)     anova(m)      Results of the code above [ **updated** ]: **anova.ez** $ANOVA                          Effect DFn DFd         F          p p<.05          ges     2           group   1  18 2.6854464 0.11862957       0.1294475137     3       direction   1  18 0.9160571 0.35119193       0.0001690471     4 group:direction   1  18 4.9169156 0.03970473     * 0.0009066868      **lmerTest::anova(model)**               Analysis of Variance Table of type 3  with  Satterthwaite      approximation for degrees of freedom                     Df Sum Sq Mean Sq F value Denom Pr(>F)     group            1  13293   13293  2.6830    18 0.1188     direction        1   1946    1946  0.3935  5169 0.5305     group:direction  1  11563   11563  2.3321  5169 0.1268      **anova(m)**               Analysis of Variance Table          Response: rt                       Df   Sum Sq Mean Sq  F value Pr(>F)         group              1  1791568 1791568 242.3094 <2e-16 ***     direction          1      728     728   0.0985 0.7537         group:direction    1    12024   12024   1.6262 0.2023         Residuals       5187 38351225    7394                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Is there a way to disable the parameter tuning (grid) feature in CARET?

text: CARET will automatically use a pre-specified tuning grid to build various models before selecting a final model, and then training the final model on the full training data. I can supply my own tuning grid with only one combination of parameters. However even in this case, CARET "selects" the best model among the tuning parameters (even though there is only one in this case), and then fits the final model to all the training data. This is an extra step I'd like to avoid. How do I simply skip the model search step across variations in the tuning grid and force CARET to build on all the training data (other than calling the underlying model library directly)?